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SUMMARY 


New methods are presented for the analysis of transient heat 
flow in complex structures, leading to drastic simplifications in 
the calculation and the possibility of including nonlinear and sur- 
face effects. These methods are in part a direct application of 
some general variational principles developed earlier for linear 
thermodynamics.'~* They are further developed in the par- 
ticular case of purely thermal problems to include surface and 
boundary-layer heat transfer, nonlinear systems with temper- 
The concepts of 


ature-dependent parameters, and radiation. 
thermal potential, dissipation function, and generalized thermal 
force are introduced, leading to ordinary differential equations 
of the Lagrangian type for the thermal flow field. Because of the 
particular nature of heat flow phenomena, compared with dynam- 
ics, suitable procedures must be developed in order to formulate 
each problem in the simplest way. This is done by treating a 
number of examples. The concepts of penetration depth and 
transit time are introduced and discussed in connection with one- 
dimensional flow. Application of the general method to the 
heating of a slab, with temperature-dependent heat capacity, 
shows a substantial difference between the heating and cooling 
processes. An example of heat flow analysis of a supersonic wing 
structure by the present method is also given and requires only 
extremely simple calculations. The results are found to be in 
good agreement with those obtained by the classical and much 


more elaborate procedures. 


(1) INTRODUCTION 


hn ADVENT of supersonic flight has brought a new 
importance to heat conduction problems in engi- 
neering. The temperature ranges involved and the 
highly transient character of the phenomena require 
anew frame of reference and the development of 
methods more suited to the new problems. 

It is our purpose here to initiate the development of 
anew approach to heat flow analysis. The expression 
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heat flow is used in its broad sense and encompasses 
heat conduction and heat transfer, convection phe- 
nomena, and radiation. The substance of the method 
has been outlined in earlier publications on irreversible 
thermodynamics by this writer.'* The principles 
were given more specific formulation in reference 3 in 
the particular field of thermoelasticity and heat con- 
duction. In those papers general methods are out- 
lined by which the elasticity and the thermal conduc- 
tion problem are treated in a unified way. The thermo- 
elastic response to thermal excitation is considered to 
result from the application of generalized thermal forces, 
defined in the same way as the mechanical forces and 
leading to Lagrangian equations for the coupled elastic 
and thermal coordinates. The thermostatic and ther- 
modynamic properties are completely defined by a 
thermoelastic potential and a generalized dissipation 
function. 

The scope of the present study includes first a re- 
formulation, in the special domain of thermal flow, of 
the previously established principles and methods. 
We make use of the concept of thermal force and of the 
dissipation function defined as previously except for a 
constant factor. We also introduce a thermal poten- 
tial. We further extend these principles to include 
boundary-layer conduction, nonlinear phenomena with 
temperature dependent parameters, and radiation. 
In addition, several entirely different methods of apply- 
ing the general principles are presented, by the treat- 
ment of specific examples, in an attempt to develop 
the art of solving complex problems of thermal flow 
analysis by procedures best suitable for each type 
of problem. It is, of course, not possible in such a 
short space to examine the innumerable variations and 
refinements of the method. We have therefore picked 
a limited number of examples and illustrations of a 
typical treatment. The heat flow in extremely com- 
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plex structures is open to analysis by such methods 
without undue complications and with suitable flexi- 
bility with regard to the introduction of intuitive and 
experimental knowledge into the formulation. 

While the mathematical basis of this new formula- 
tion is presented in the variational language, it is not 
a straightforward application of the variational calculus 
to the heat conduction equations. The physical phe- 
nomenon is considered to be represented by a vector 
field instead of a temperature field. This opens the 
way to new concepts such as generalized thermal force 
and to straightforward procedures for the interconnec- 
tion of thermal networks. 

From the standpoint of simplicity, compared to the 
classical treatment of heat conduction, the new ap- 
proach bears similarity to the powerful simplification 
which was achieved in the vibration analysis of com- 
plex structures by the introduction of Lagrangian 
methods. We must bear in mind, however, that the 
steps are far from identical because of the essential 
physical difference between the nature of the 
ena of diffusion and dynamics. One of our 
here is to show how special adaptation of the methods 
to heat flow phenomena lead to remarkable simplifi- 
cations. An unexpected result is a simplification which 
leads to a nonlinear differential equation with a single 
Another 


phenom- 
purposes 


unknown even in a physically linear problem. 
special feature of interest lies in the use of time- 
dependent approximate solutions, instead of the usual 
static distributions which are of standard use in the 
Rayleigh-Ritz solutions of dynamical problems. This 
also leads, of course, to a method of successive approxi- 
mations, whereby a previous solution may be intro- 
duced in the following step and ordinary differential 
equations established for the correction. This con- 
stitutes, by the same token, a procedure for testing the 
accuracy of an approximate solution. 

The basic principles in variational form as applied to 
a thermal flow field are introduced in Section (2). 
They are a direct consequence of the more general 
theories developed in references 1, 2, and 3.  Itis shown 
how boundary-layer heat transfer may be included. 

The description of the flow field by means of general- 
ized coordinates and the concept of thermal force are 
introduced in Section (3). By means of a thermal po- 
tential and a dissipation function, it is shown how the 
thermal flow problem may be formulated by differential 
equations of the Lagrangian type. The inclusion of 
surface and boundary-layer heat transfer is also dis- 
cussed. Section (4) deals with linear systems and dis- 
cusses thermal modes, normal coordinates, orthogo- 
nality, and general expressions for the thermal ad- 
mittance in operational form. The basic mathematical 
and thermodynamic background for these concepts and 
their properties was established in references | and 2. 
Various ways of applying the general equations to ther- 
mal flow problems are discussed in a general vein in 
Section (5). Specific procedures for one-dimensional 
problems are developed in Section (6). One example 
uses normal coordinates and a modification of the latter 


AERONAUTICAL 


1957 


DECEMBER, 


SCIENCES 


for flows deviating only slightly from the steady stg, | 
Simpson's method of integration is used to write it 
equation directly in matrix form particularly syjtul 
to solution by methods of iteration, as in vibratin} 
analysis. A particularly simple way of analyzing he, 
conduction in a slab is illustrated in Section (7). Thy 
leads to a simple differential equation with one y 
known. Comparison of the results with the exay 
series solution shows the method to be warps 
accurate considering its simplicity. ’ 

The principles are further extended in Section (S) 
include systems where dependence of the heat capacit 
and conductivity on the temperature leads to nonlines 
equations. Temperature-dependent heat-transfer « 
efficients as well as surface radiation conditions ar 
included in the formulation. This is illustrated j 
Section (9) by the example of heat conduction in a sla 
with temperature-dependent heat capacity. It 
found that, because of the nonlinearity, the coolin: 
proceeds quite differently from the heating. 

As a final example, in Section (10) we treat the prob 
lem of a typical portion of a supersonic wing structur 
under aerodynamic heating. The example was solve | 
previously by Pohle and Oliver.6 The present results| 
obtained with very little effort, are compared wit! 
those calculated by the exact and considerably mor, 
elaborate procedure, and their accuracy is found to by 
Various methods of refinement © 


quite satisfactory. 
the solution of this particular problem are also dis 


cussed. 


(2) FUNDAMENTAL VARIATIONAL PRINCIPLES IN HE. 
CONDUCTION 


Consider an isotropic body, with a thermal condue- 
tivity k(x, vy, s) and a heat capacity c(x, y, 2) per unt 
volume as functions of the coordinates. For our pur 
pose it is convenient to introduce an excess temper 
ature 6 = 7 — 7) over an equilibrium temperature / 
The temperature @ satisfies the equation of heat con- 
duction, 

(0 Ox) [k(08 Ox)}] + (0/Oy) [R(08/Ov)] + 
(0/0z) [k(00/0z)] = Of) (2.1 

We shall formulate a variational principle which 
equivalent to this equation. To this effect we intro- 
duce a heat flow vector field H such that the rate « 
heat flow at every point is OH/Of per unit area, norm 
to H. Energy conservation is expressed by the relatio! 


co = —div H (2.2 
We define a thermal potential 
V = (1/2) fff (2.3 
obtained by integration in the volume 7. We als 
define a variational invariant 
(2.4 


6D = k) (OH/dt) -6Hdr 
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We shall see that |’ plays a role analogous to a poten- 
‘tial energy while 6D is related to the concept of dissipa- 
tion function. Except for a constant factor they are 
identical with the functions introduced by the writer 
in references 1, 2, and 3 in connection with linear irre- 
versible thermodynamics. It was shown in these 
-eerences that I’ may be considered as a generalization 
of the free energy for systems which are not in thermal 
equilibrium. The relation of 6D to a dissipation func- 
tion will appear more clearly in the next section dealing 
with generalized coordinates. 

With these definitions we state the variational prin- 


6V + 6D = ff én-6HdS (2.5) 


| iple as 


This equation is to be verified identically for all 
arbitrary variations 6H of the heat flow field in the 
volume 7. The surface integral is extended to the 
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boundary S of the domain with the unit normal n 
to the boundary taken positive inward. To show that 
Eq. (2.5) is equivalent to the heat conduction equation, 


we evaluate 


= 66007 = 66(div H)dr (2.6) 


Integrating by parts, we find 


v= ff | 6H- grad + 6n-6HdS (2.7) 


Substituting Eqs. (2.4) and (2.7) in Eq. (2.5) and 
putting the coefticient of 6H equal to zero yields 


grad 6 + (1,/k) (OH, Ot) = O (2.8) 
From Eqs. (2.8) and (2.2) 

div (k grad 6) = c(06/Ot) (2.9) 
which is identical with Eq. (2.1). The validity of the 
variational principle (2.5) is thus established. 

The above derivation assumes a body of isotropic 
thermal conductivity. In the more general case of 
anisotropy the thermal conductivity is defined by a 


symmetric tensor 


= (2.10) 
and the law of thermal conduction is written 
J 
k;,(00/0x;) = —(OlI,/dt) (2.11) 


Eq. (2.1) is replaced by 


(0/dx,) [k,,(00/dx,)] = c(00/dt) (2.12) 


Proceeding as above for the isotropic case, it can be 
shown that the variational principle (2.5) is also valid 
lor the general anisotropic case, provided 6D is replaced 
by 


ij 
5D = (OH, /0t)-6Hjdr (2.13) 
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In this expression \,,; is the inverse matrix of k,; and 
represents the thermal resistivity 


Ay = [Ri] (2.14) 


The law of thermal conduction is then 


00/d0x, = —> (2.15) 


The principles as stated above may be readily ex- 
tended to include the effect of surface and boundary- 
layer heat transfer. This can be seen immediately 
since the heat-transfer layer is equivalent to a material 
with zero heat capacity. This will be developed more 
explicitly in the next section. The question of radiation 
loss will be discussed in Section (S) in connection with 
nonlinear problems. 

From an intuitive viewpoint, it is interesting to point 
out the complete analogy of heat conduction and the 
seepage of a compressible viscous fluid through a porous 
solid. Such an analogy was discussed in reference 3. 
The mass flow rate corresponds to the rate of heat 
flow 0H O/, the pressure to the temperature and the 
increase of fluid mass per unit volume to the heat 
content. The fluid compressibility represents the 
heat capacity and the permeability is the equivalent 
of the thermal conductivity. The analogy holds for 
nonisotropic media and may be extended to nonlinear 
systems, as will be demonstrated in Section (8). It 
also indicates that such concepts of mechanics as gener- 
alized forces and coordinates may be introduced in 
thermal problems. It can be done directly without 
referring to the mechanical model, as will be shown in 
the following section by applying the variational prin- 


ciple demonstrated above. 


(3) GENERALIZED COORDINATES AND THERMAL 
FORCE—-SURFACE HEAT TRANSFER 


We now introduce the concept of generalized coor- 
dinates in conjunction with the variational principle. 
There are many ways in which these generalized coor- 
dinates may be defined. We may express the flow 
field H as a sum of fixed configurations H,(x, y, 2) with 
a variable amplitude of each configuration 


H = >) (3.1) 


The unknown amplitudes g,; are the generalized coor- 
dinates. However, these coordinates need not be re- 
lated linearly to the field H. More generally, we could 
write 


H= 


In this case the generalized coordinates are simply a set 
of m parameters defining the field configuration. As 
we shall see in the examples [Section (7)| it is advan- 
tageous sometimes even in linear problems to choose 
generalized coordinates which are related nonlinearly 
to the unknown thermal variables. 

Consider now the variational relation (2.5) and let 
us apply it first to an isotropic medium. 


7 
4 
(2.3 


variations of g;, 


6H = > (0H/dq,)6q; (3.3) 
Also we may write 
(0/o))H = >> (0H/0q))4; (3.4) 
Hence, 
(0/0q,;) (OH /Ot) = 0H/0q; (3.5) 


The invariant (OH/Ot)-6H in Eq. (2.4) may, therefore, 
be written 


(0H/dt)-6H = (0H /dt)- >> (0H/0q,)5g; = 


(0H/dt) (0/dg,) (OH /dt)ég; (3.6) 
or 


(0H/dt)-6H = 6q,(0/dg,) [(1/2) (OH/dt)?]_ (3.7) 


With these results the variational equation (2.5) be- 
comes 


+ = f f on-6HdS (3.8) 


with 


0D /04g; = ‘k) X 


(OH/Ot)-(0/0g;) (OH/Odt)dr (3.9) 


Hence, we may introduce the invariant 
D = (1/2) (1/k) (OH/dt)? dr (3.10) 


Finally, we put 

6n-d6HdS = (3.11) 
with 

0; = ff (OH /0g;)dS (3.12) 


6V = (3.13) 


The variational principle, therefore, yields n equa- 
tions for the field parameters g;—namely, 

(OV/0g:) + (OD = Q: (3.14) 

These are the analog of the Lagrangian equations for 

a mechanical dissipative system with a potential energy 

V and a dissipation function D. The above equations 

have been derived for an isotropic medium. In the 

anisotropic case, it is easily verified that the same 
equations hold provided we define D as 


D = (1/2) f f (OH,/dt)dr (3.15) 


This expression for D as well as Eq. (3.10) shows its 


860 JOURNAL OF THE-AERONAUTICAL SCIENCES—DECEMBER, 1957 


The variations 6H in this case are due entirely to the 


physical significance as it is related to the rate of » 
tropy production.’ We shall refer to D as a dissip, 
tion function also in purely thermal problems. 

The generalized force Q; will be referred to as th, 
thermal force. It can be seen that Eq. (3.11) defing 
it in exactly the same way as a mechanical force—j¢ 
as the work done by a temperature @ on a virtual dj 
placement 6H. 

The above derivation does not require that the fieli 
H be completely defined by the coordinates 4g, as ;; 
Eq. (3.2). It is possible to have the time appex 
explicitly in the description and write 


= .. . du; ¥, %, (3.11 


This will lead to the same differential equations (3.14 
for the coordinates g; but with the time variable appexr. 
ing explicitly in the equations. This point is of in. 
portance in applications as it leads very often to cop. 
siderable simplification to adopt a_ time-dependent 
description of the field. 

Let us now consider the case of a surface heat. 
transfer effect and show how it may be included in th 
above formulation. We denote the temperature out. 
side of the layer by 6, and the temperature of the sur. 
face of the body by @. With a heat-transfer coef. 
cient A for the layer, the local rate of heat flow acros 
the layer is 


OH,,/ot = — @) (3.17 


The normal component of the heat flow into the bod; 
is denoted by //,. From Eqs. (3.12) and (3.17) the 
thermal force at the surface of the body is 


Q; = ff = (1/K) (O/7,,/dt) X 


+ ff (3.18 
Because of Eq. (3.5) we may write 


oD, /d4; = f (1/K) (0H,,/dt) (OH, /dq,)dS = 


ff (1/K) (O77,,/Ot) (0/0g;) (OH, /Ot)dS (3.1! 


With a surface dissipation function, 


D, = (1/2) ‘K) (O0H,,/ot)*dS (3.20 


Hence, Eqs. (3.14) may again be applied to the tote 
system including the surface layer provided we inclut: 
D, in the total dissipation function and write 


Di =-(1/2) fff (1/k) (OH/ot)*?dr + 
» [fo ‘K) (OH,/ot)*dS (3.2! 


We must also use a definition of the thermal force l) 
means of the temperature outside the layer and wnt 
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0, = 6, (OH, /dq,)dS (3.22) 


When the surface layer is a moving fluid boundary 
aver, , is the so-called “adiabatic wall temperature.” 
It is of interest to point out that in Eq. (3.21) the 
poundary heat-transfer coefficient A may be time de- 
pendent. This is particularly useful in problems of 
heat transfer through boundary layers with variable 


fuid velocity. 
(4) THERMAL ADMITTANCE AND NORMAL COORDINATES 


The case where Eqs. (3.14) are linear is of particular 
interest. This will happen, of course, if the physical 
problem 1s linear and if the generalized coordinates are 
related to the thermal field by a linear expression of 
the type (3.1). The thermal potential V’ and the dis- 
sipation function D are then positive definite quadratic 


forms 
V = (1/2) > a9 
(4.1) 
D = (1/2) 


and Eqs. (3.14) read 


j 
aug) + bugs = Q: (4.2) 
with symmetric matrices, = @);, = 
We consider first the case where no thermal forces are 
applied (0; = 0). Solution of the homogeneous equa- 


(3.17) the|tions requires the solution of the characteristic deter- 


IT, X 


dS (3.18 


S = 

1S (3.19 
S (820 


minant 


det\a,, + pb,,| = 0 (4.3) 


with p as unknown. It can be shown that the roots 


~\, are real and never positive.'! The characteristic 


solutions are 


= ™ (4.4) 


where the ¢,‘°’ are normalized modal distributions and 
C” arbitrary constants. We refer to these character- 


stic solutions as thermal modes. These modes satisfy 
thogonality relations (r s) 
ij 
a)” = = 0 (4.5) 
lhe proof of these relations is identical with that given 
nreference | in connection with relaxation phenomena. 
Whether the roots \, are zero or infinite or multiple, 


» the tote) there are always n such orthogonal modes in a system 


we include} , 


fn coordinates. 
These mathematical properties are identical with 


those of a dynamic system with inertia and elasticity. 


The well-known numerical methods used in vibration 
alysis, such as iteration, for instance, may be applied 


1S (3.21 to the evaluation of the thermal modes. 


1 force bi 
and write 


We should note that, in general, a thermal system 


will contain heat flow configurations with conservative 
ields—i.e., such that div H = 0. The thermal po- 
tential |” is independent of the particular generalized 
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coordinates corresponding to these fields, while D de- 
pends on these coordinates. These fields also consti- 
tute thermal modes of infinite relaxation time—i.e., of 
zero roots\,. The case is analogous to that of ignorable 
coordinates in dynamics, such as, for instance, the 
coordinate corresponding to the free motion of a solid. 

Once the thermal modes are known, the basic equa- 
tions (4.2) may be expressed in terms of normal coordi- 
nates defined by the relations 


a= (4.6) 


With these coordinates all equations are uncoupled 
since |” and D reduce to sums of squares. Using 
properly normalized modal distributions, Eqs. (4.2) 
become 


dE, + & = E, (4.7) 


The generalized forces are 


z= 9,0, (4.8) 
The operational solution of Eq. (4.7) reads 


= =/(p + d,) (4.9) 
with p = d/dt 


Substituting Eqs. (4.8) and (4.9) into Eq. (4.6) we 
obtain the operational solution of the original Eqs. 
(4.2) as 


a= LAR, (4.10) 


with an operational matrix 


Ay =Ay = [Cy/(p + d,)] (4.11) 


which constitutes the thermal admittance of the system. 


The coefficients are 


Cy = (4.12) 


If there are multiple roots, we may collect all the terms 
corresponding to the same root into one single term 
so that the summation in Eq. (4.11) for C,* may be 
taken to extend to all distinct roots only. 

In a thermal system there are no infinite values of 
\,—i.e., all relaxation constants are finite. However, 
as pointed out above there may be zero roots \, corre- 
sponding to steady heat flow through the system. This 
may be expressed by separating the zero root term in 
the admittance, 


Ay = [Cy/(p + AD) + (Co/p) (4.13) 


The significance of these operators is brought out by 
considering the case where all thermal forces vanish 
except Q;. Assuming the time dependence such that 
Q; corresponds to a sudden rise in temperature at / = 0 
expressed by a Heaviside step function, we write 

(0 t<0 


QO; = = (4.14) 
(1 


_ ; 
| 


862 JOURNAL OF THE AERONAUTICAL SCIENCES DECEMBER, 1957 
The coordinate gq; is then The subscript » indicates normal components of th, out, th 
C0) C, vectors. Now we may choose the fields F; and @, represt 
+ 1(t) (4.15) such a way that 
From the significance of these operators! we obtain bi’ = fff + 
s of unk 
= (Cy/ds) + Cot (4.16) (1/K)FOjndS = 0 (55 eral 8) 
The bracketed term corresponds to the establishment ierent 
of an equilibrium temperature distribution, while the i.e., such that the @,'s are orthogonal to the divergeng.} norm 
term proportional to / represents the steady-state heat free field F;. The thermal potential I” is independey — 
flow through the system after the equilibrium temper- of We have to fiel 
atures are established. 
Vv = (1/2) fff (1/c) (div H)*dr (5.9 
(5) Some GENERAL PROPERTIES AND PROCEDURES . The v 
There are certain general properties of the temper- = (1/2) = 
atures and thermal flow fields which are worthy of a 
closer scrutiny because of their bearing on basic meth- with dj; “fff. (1/c) (div @;) (div @))dr (5.1 The 
ods of solution of the equations. and ¢ 
It was pointed out above that, in general, a thermal Under these conditions the equations for g; are inde. 
system will contain a divergence-free field of flow cor- pendent of the coordinates /,, and we may write 0, = 
responding to ignorable coordinates. Let us examine j j In tet 
this point more closely by separating the flow field into aig) + bugs = Qi (9.12) condi 
a divergence free field Hy and a remainder Hy, The question arises of how to establish such a cooréi- 
H=H,+H, (5.1) nate system. This may be done by introducing scalar 
field distributions y,; associated with a thermal flow 
The field H, is associated with temperature changes field by This 
while H, is not. Each field may be represented by its _ | ortho; 
own set of generalized coordinates g; and /; = k grad (9.13)" obtaii 
Then by Green's formula, and taking into account} equal 


He = Om: (5.2) 


The /;’s represent the ignorable coordinates. By defi- 


nition, 
div H, = div F; = (5.3) 
Since the temperature is 
= —(1/c) div H = —(1/c) div (5.4) 
the thermal potential |” is independent of the /; co- 


ordinates. The dissipation function including the sur- 


face heat transfer is 


or = (1/2) by + 
(1/2) dy + (1/2) > (5.6) 
with 
1 1 
bi; = [ff ff 
bis = fff (5.7) 
1 1 
be =f povodrt J 


Eq. (5.3), we find 


where is the normal to the boundary taken positive 
Condition (5.8) becomes 


Fin[W; + (1/K)OjnJdS = 0 (5.15 


This is satisfied if the scalar y; satisfies the following 
condition at the boundary: 


+ (k/K) grad, = 0 


Because of Eq. (5.3) we may also substitute a constant 
on the right-hand side of this condition instead of 
zero. With such a coordinate system we have to solve 
the differential equations (5.12) for g;, each of which is 
associated with a particular temperature configuration. 
The equations for /; are irrelevant as far as tempera- 
ture is concerned and are only of importance if we wish 
to calculate the heat flow. The equations for /; are 


(5.14 


outward. 


(5.16 


> = (5.17) 


where Q,” is the thermal force associated with /;. The 
solution of this system is trivial. As we have pointed 


* If kis discontinuous, we must choose y¥; in such a way that the 
normal component k grad, yi is continuous across the surface of 
discontinuity of k, as well as the corresponding temperature. 
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out, these ignorable coordinates are the analog of those 
representing the free body motion in the vibration anal- 
ysis of an elastic structure. They are eliminated here 
py a similar condition of orthogonality. 

“The differential system (5.12) with a reduced number 
of unknowns possesses the same properties as the gen- 
eral system with its own characteristic values, all dif- 
jerent from zero. If the system is represented by its 
normal coordinates, the matrices and are diago- 
nalized. Two normal coordinates £, and £, correspond 
to fields 


H,” = @%:,, = @%, (5.18) 


The vanishing of the nondiagonal term a,, means 


Ay, -| c) (div 8”) (div = 0 


The temperature fields associated with the modes £, 


(5.19) 


and &, are 


= —(1/c) div 0™, @, = —(1/c) div (5.20) 


In terms of the temperature fields the orthogonality 
condition (5.19) becomes 


{ff 0dr = 0 (5.21) 


This condition is valid for nonisotropic media. The 
orthogonality condition in terms of the flow fields is 
obtained by putting the nondiagonal coefficients 6,, 
equal to zero. We find 


+ 
{f (1/K)0,0,dS = (5.22) 


For an anisotropic medium it is easily verified by using 
Eq. (3.15) that the orthogonality condition becomes 


[| + 


Sf (1/K)0,”0,@dS = 0 (5.23) 


where 0,” and 0,“ are the Cartesian components of 
the flow field. 

In solving Eqs. (5.12) for the coordinates we may 
first calculate the normal coordinates. The normal 
coordinates may be conveniently obtained by an iter- 
ation procedure applied to the homogeneous system 
Le, putting Q; = 0. The homogeneous equations in 
matrix form are 


lay] = [bus] (5.24) 


For an exponential solution of the type e 
(1/d) [gi] = Las) 


The iteration procedure is the same as the well-known 
one used in vibration analysis. Substituting a column 
4; on the right-hand side, we obtain a new column 
Which is again substituted. This process converges 


we write 


(5.25) 


toward the thermal mode of smallest value of the re- 
laxation constant ». Use of the orthogonality condi- 
tion yields successively the modes of higher i. 

We notice that if we use the thermal modes as the 
coordinates, the steady-state solution is represented 
by a superposition of a complete set of modal distribu- 
tions. 

In many problems with slowly varying temperatures 
the deviation from the steady state may be small and 
only a small number of relaxation modes with slow 
decay are excited. In such a case it is preferable to 
Assume we have calculated the 
The equations 


proceed as follows. 
normal coordinates of the system. 
then are 


rt, + & = =, (5.26) 
If the forces were applied very slowly, the solution for 
£, would be 


(5.27) 


> 


We may separate the solution into an instantaneous 
steady state ¢, and a correction £,—1.e., replace &, by 


+ 


This yields 


(5.28) 


+ & = —& = 


The force =, is found by calculating the virtual work of 
the time derivative of the boundary temperature 6 on 
each relaxation mode. It is seen that if the tempera- 
ture varies slowly, only a small number of modes are 
excited. This is particularly advantageous if we evalu- 
ate the modes by iteration since only a few modes of 
lowest value of \, need be calculated. The steady-state 
solution ¢, is not evaluated by normal coordinates, but 
the temperature field may be calculated directly. 
This is best illustrated by the example in the following 
section. 

Before closing this discussion, attention should be 
called to an important feature. With the exception of 
normal coordinates the methods discussed above are 
not restricted to component field configurations ©, 
which are constant. We may introduce variable fields 
such as 


1 
H, = @(x, y, 2,0) + Ox, 2, (5.29) 


This is extremely important since it opens the way to a 
method of successive approximation where @» is an 
approximate solution to which qg, terms are added as a 
correction. In this case the differential equations for 
g; may or may not have time-dependent coefficients. 

Another important point, already mentioned before, 
is the possibility of including boundary heat-transfer 
coefficients A which are time dependent. Such is the 
case, for instance, in aerodynamic heating. 

Methods outlined here all lead to linear equations. 
They are not the only ones available. In some prob- 


lems it is preferable to choose the generalized coordi- 
nates in a way different from the above. This will be 
illustrated in Section (7). 


= 0) (5.4 
(5,9 
| 
(5.15 ; 
td 
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(6) APPLICATION TO SOME ONE-DIMENSIONAL 
PROBLEMS 


We shall now illustrate the methods of thermal flow 
analysis, outlined in the previous section, by formu- 
lating them more precisely on some problems of one- 
dimensional flow. Other methods relative to one- 
dimensional flow will be introduced in the following 
section. 

In the present approach the unknown variables rep- 
resent the amplitudes of certain given field configu- 
rations. In reference 3 we have discussed the solution 
of a plate problem by means of normal coordinates. 
This was the problem of a plate brought suddenly to a 
temperature #) on one face and kept at temperature 
6 = (0 on the other. Here we shall discuss the case 
where the plate is thermally insulated on the other 
face. The temperature is represented as a Fourier 
series, 


h = = >> q, sin [n + (1/2)] 
) 


( 


(6:1)* 


and the corresponding heat flow is derived by integra- 


tion. Since // = Oat y = /, we write 


cos 

0 + (1/2)] 27 

This corresponds to distributions of the scalar y,, 


= sin [n + (1/2)] (ry//) 


(6.2) 


(6.3) 


satisfying condition (5.16). 
The thermal potential and dissipation function are 


2c Jo 4c 0 | 


oH 2 [3 Gn? 
= | 


These forms reduce to sums of squares because the 
component field configurations satisfy the orthogonality 
relations (5.19) and (5.22). The thermal force Q, 
-ach coordinate is derived from the 


associated with 
variation 6/7 due to 6g, at y = 0, 


= 061 = (6.5) 


a(n + (1/2)] 
Hence, the differential equation (3.14) in the present 
case takes the form 


l 


These are uncoupled equations as should be since we 
have chosen normal coordinates. The equations are 
valid, of course, for any arbitrary variation of 0) with 


time. In our particular example, we have 


* These sine functions are chosen because they are the natural 
thermal modes of the system. It can be inferred from results 
of the previous section that they constitute a complete orthogonal 


set of functions. 
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= (6.7) 
where 1(f) is the Heaviside step function, 
(jo 
I(t) = 
(1 ¢>0 


In more general cases the choice of normal coordi. 


nates as above will have the inconvenience of the uni. } 


formly distributed temperature at equilibrium being 
represented by a slowly convergent Fourier series. |} 
may be preferable therefore to proceed as already sug. 
gested in the general discussion of Section (5). The 
temperature is represented by a term 6(¢) independent 
of y, corresponding to the instantaneous steady state 
and a residue (;, 


0 = A(t) + A (6.8 


We apply the procedure leading to Eq. (5.28). The 
residue 6, is represented by normal coordinates 4,, satis- 
fying the same Eq. (6.6) except that the generalized 
force is now due to 4) and is divided by the character. 
istic value 


An = (Rw?/cl?) [n + (1/2)]? (6.9 
The equations are 


/ / | 
Qn + qn => Ay 


(6.10) 
2 22, + (1/2)] 


In the present case, where % is given by Eq. (6.7), the 
solution is 

9 

2c 


~ e (6.11) 
+ (1/2)] 


The temperature 6; associated with this coordinate is 
given by Eq. (6.1) 


2c 1\ ry 
e sin (» ) — (6.12) 
0 min + (1/2)] l 
The total temperature field is 
+0 = — > a(n 2)}} x 
0 
sin [n + (1/2)] (ry/l) (6.13) 


which is the classical solution for this problem. 

It is, of course, possible to use any type of function, 
such as trigonometric functions, which are not neces- 
sarily orthogonal, and polynomials. For instance, we 
could use a polynomial 


x 


0 


(6.14) 


This might be particularly useful for a heterogeneous 
material, such as a slab made up of layers of different 
materials. 

A convenient method for numerical work is to take 
advantage of the accuracy of Simpson's rule in pet- 
forming the integrations by adapting it to the present 
The slab is divided into an even number, 2, 


case. 
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of intervals and the temperatures at the point of sub- heat flow is obtained by integration of the parabolic 
In pairs of intervals area.* By integrating‘ in the successive intervals, the 
heat flow vectors //; at the various points of sub- 


(6.7) 
division are = O, 1,...2N). 


| -5 12-23 1/12 0 0 
lal coordi. H. — H, | 112 —2/3 —5/12 0 0) 6; 
the uni. i=; 0 1/12 (6.15) 
um being cay | 0 1/12 —2/3 —5/12 
eady sug. 
(5). The & by addition of successive lines 
I, — —5/12 -2/8 1/12 0 0 
| Ho -1/3 -4/8 -1/3 0 0 
=; —1/3 —4/3 -—3/4 -—2/3 1/12 (6.16) 
(68 -13 -4/3 -23 -43 -13 
Boy — Hy 
S). The 
Satis- By putting = cAy 
neralized 
" this may also be written f 
aracter- 
q-1 
4, 
(6.9 1 ‘ 
. (6.17) 
cAy 
Hoy Boy 
6.10 
0 0 0 0 = 
7). th 1 —8/12 —2/3 1/12 0 0 
the 9 9 > 
1 —1/3 —4/3 -—1/3 0 0 
itl 6.18 
[4] 1 -1/3 -4/3 -3/4 -2/3 1/12 — 
1 —1/3 —4/3 -4/3 -13 — 
(6.11) | 
nate js | [he dissipation function is evaluated by Simpson's method expressed in matrix form 
. . 
D = (12k) (O77 Ot)*dy = (Ay/2k) [S| | (6.19) 
(6.12) 
where (Wi) = 
= (Ho... (6.20) 
ction, and [S| is the diagonal matrix corresponding to Simpson’s rule 
e, we i 4 0 
1 4 
6.14) [S] = = (6.21) 
3 
eous 0 
‘rent 4 
take 
per- 
sent * This integration of a continuous distribution of temperature in order to calculate the coefficients of the differential equation is the 
2N, reason for the considerably increased accuracy of the present method over the usual finite difference methods of solution. 


the temperature is represented by a parabola passing 


) through the three ordinates of the two intervals. The 


division are given by 


865 
| 


JOURNAL OF 


THE AERONAUTICAL SCIENCES—DECEMBER, 


VAN 
\\\ 


\\\\ 


h 
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Distribution of h = Cé for first phase (1), and second 
phase (2), in the heating of a slab. 


Fic. 1. 


Substituting in the expression for D 
D = (c?/2k) (Ay)* [B] [4:] 
| 


(6.22) 


with the matrix = [g9:] (6.23) 


bien! 
We have introduced the transposed [g;]’ of [g,] and 
[B] = [A]’ [S] [A] (6.24) 


where [A |’ is the transposed of [A]. 
tential is again by Simpson’s rule 


The thermal po- 


V = (c/2) = (c/2)Ay[6;]’ (6.25) 
We may consider the 2V + 2 parameters g-1, 6... 
62x as the generalized coordinates and proceed as in the 
general case. We must evaluate the corresponding 
generalized forces. Let us assume that the temper- 
ature #2, is given while the surface at 7 = 0 is insulated. 
This means Hy) = 0. In this case there are 2V + 1 
coordinates 69; . . . Oey and the matrix [A] is reduced 
to that in relation (6.16). The generalized forces are 
evaluated by the expressions 


Hon = QO 66; (6.26) 
The differential equations are 
[S] [0:] + (c/k) (Ay)? [B] [4:] = [Q:] (6.27) 


Since [.S] is diagonal, the therinal modes may be ob- 
tained immediately by iteration. 

This method is quite general. If we are dealing with 
a composite slab, the procedure is identical. 

The case where there is a coefficient of heat transfer 
K at the surface of the slab is easily taken care of by 
adding a term to the dissipation function in accordance 
with Eq. (3.21). 

The additional term is 


D = (1/2K) Hey* 


The thermal potential remains unchanged. 


(6.28) 
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(7) HEATING OF A SLAB WITH CONSTANT 
PARAMETERS—THE CONCEPTS OF PENETRATION 
DEPTH AND TRANSIT TIME 


In the previous section we have outlined procedures 
for the analysis of one-dimensional problems leading 
linear differential equations for a sequence of unknowy 
parameters representing the amplitudes of component 
flow fields. The flexibility of the present method is wel] 
illustrated by the possibility of using an entirely differ. 
ent approach where only one unknown parameter jg 
used to represent the unknown field. Although the 
problem is physically linear, the parameter is found to 
satisfy a nonlinear differential equation of the first 
order. The nonlinearity in this case is compensated 
by the simplicity of the field representation. This 
approach also leads to two very useful concepts in 
practical applications—1i.e., penetration depth and transit 
time. 

Consider a plate of thickness / with constant values 
of the thermal conductivity k and heat capacity c (see 
Fig. 1). 

One face at y = O is heated suddenly to the tem- 
perature 6) at ¢ = 0. The other face at y = / is ther- 
mally insulated so that no heat flows across it. 

We shall assume the temperature distribution to be 
represented by a parabola and consider two phases in 
the phenomenon. In the first phase, the temperature 
has not yet begun to rise at the opposite wall and 
everything occurs as if the wall thickness were infinite. 
During this phase the heat content distribution 4, 
which is proportional to the temperature, is approxi- 
mated by 


h 


ho [1 — for vy¥<uq 
(7.1) 


for 
hk = cb, 


yon 


with ho = c% 


The parameter 4; is the generalized coordinate and may 
be called the penetration depth. The problem is to es- 
tablish the value of g; as a function of time. 

The thermal potential is 


gi 
» dy = (1 29 h*dy (7.2) 
0 0 


hence, V = (1/10) (Ao?/o)q (7.3) 


In order to determine the dissipation function we 
must introduce the heat flow vector H/. This is ob- 
tained from by integrating (2.2)—i.e., from 


—dH/dy =h (7.4) 


Taking into account the condition J/7 = 0 at y = q; we 
find 


H = (1/3)hog: — hoy + ho(y?/q:) — ) 


(ho/3) 
> (7.5) 
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The dissipation function ts 


D = (1 2k) [ (H)*dy = 
0 
(1/2)-(13/315) - (he? (7.6) 


Finally, we must evaluate the generalized force. This 
is obtained by considering the virtual displacement 
jH at y = 0. From Eq. (7.5) we have 


= (1/3)hodqr (7.7) 

By definition, 
= 06H = (1/3) (7.8) 
Hence, = (1/3)ho% = (7.9) 


The differential equation for q, is 
(OV/dq:) + (OD/d4:) = (7.10) 
or (1,10) (Ao? /c) + (13/315) (Ae? = 
(1/3)ho? (7.11) 


This is a first-order differential equation for g,, which is 
easily integrated. With the initial condition g; = 0 


fort = 0 we find 
= 3.36V (k/c)t (7.12) 


In the first phase the penetration depth is proportional 
to Yt. It will be noted that we have treated a linear 
problem by a nonlinear differential equation. The non- 
linearity in this case is of a purely geometrical nature 
and is introduced by the fact that the generalized co- 
ordinate is related nonlinearly to the physical variable 
to be determined. The general disadvantage of non- 
linearity is, in this case, largely outweighed by the re- 
markable simplicity of the solution. 

From Eqs. (7.1) and (7.12) we derive the temperature 
distribution during the first phase 


y 
6 = = | (7.13) 
3.36 VR(t/c) 
This phase ends when the temperature begins to rise at 
the boundary—i.e., at a time when g, = / or 
ty = 0.0885 cl2/k (7.14) 


We shall call this the transit time. With this transit 
time we may also write 


II 


IVt/t; 
6 = [1 (y DMV t, |*\ 


In the second phase the heat content c# at y = / is 
used as generalized coordinate gz and the heat content 
distribution is again approximated by a parabola 


(7.15) 


h = (ho — qe) [1 — (v/D]? + @e (7.16) 


The thermal potential is 
V = (1/2c) h?dy = (1/2c)ho"l + 
0 


15c) (qe ho)? 3c) (q2 ho) (7.17) 
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—— PRESENT METHOD 
— —EXACT (REF. 5) 


= = 4.0 
ty 


~ 


0 
9/2 
Fic. 2. Temperature distribution at various times in a slab 
The heat flow vector is found by integrating Eq. 
(7.4) with the boundary condition /J = 0 at y = /. 
We find 


—H = (ho — qe) ly — + (y*/3P) — (1/3)) + 
gly — 1) (7.18) 
H = + (y*/3P) + (2/3)/) (7.19) 
Hence the dissipation function is 
D= (1 2%) = (34/315) (7.20) 
We must finally evaluate the generalized force Q. by 
considering the value of 6// at y = 0. We write 
= = (2/3)Al6q2 (7.21) 
Hence, Q2 = (2/3)0ol = (2/3) (ho/c)l (7.22) 
The differential equation for q» is 
(OV/Og2) + (OD /Og2) = Qe (7.23) 
or explicitly, after simplification and introduction of 
the transit time /;, 
g2 + 4.67hgGe = ho (7.24) 
The solution with initial conditions g. = 0 at f = 4, is 
(q2/ho) = 1 — exp | —0.214[(t/t)) — 1]f (7.25) 
With this value of go/h) the temperature distribution 
in the second phase is given by 
6/00 = [1 — (qe/ho)] [1 — + (7.26) 
This expression, along with Eq. (7.13), gives the com- 
plete time history of the temperature. The distribu- 


tions 6/6) thus obtained are plotted in Fig. 2 as func- 
tions of y// for three values of the time ratio 


t/t; = 0.2, 1.0, 4.0 


The dotted line shows the results obtained by the 
classical Fourier series expansion and is taken from 
reference 5. Agreement of the present method with 
the exact solution is seen to be quite satisfactory. 

It is of interest to evaluate the order of magnitude 
of the transit time /; [see Eq. (7.14)]. The value of 


6 
tr 
t=02 
t, 
0 
| 
4 
7.5) | 
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the diffusivity x = k/c is x = 0.12 cm.?/sec. for steel. 
Considering a steel plate | in. thick, / = 2.5m. and f; = 
4.6 sec. For a dural plate of same thickness x = 0.86 
cm.?/sec., 4 = 0.62 sec. 

This concept of transit time is useful in many ways. 
It gives immediately the order of magnitude of specific 
transient effects in complex structures and indicates 
which features of the thermal flow may be neglected 
in comparison with other time constants. 

The solution obtained here for the sudden rise in 
temperature at the wall involves very simple functions. 
It is, therefore, very simple to express in closed form 
the temperature field originating when the wall tem- 
perature is an arbitrary function of time using Duha- 
mel’s integral.‘ This will be applied in the example 
of Section (10). 


(S) EXTENSION TO NONLINEAR SYSTEMS WITH 
TEMPERATURE DEPENDENT PARAMETERS AND SURFACE 
RADIATION 


In a medium where the heat capacity c(x, y, 2, #) and 
thermal conductivity k(x, y, 2, 6) are functions of the 
temperature 0, the thermal conduction is governed by a 
nonlinear equation identical in form with Eq. (2.1). 
We will show that the variational principle as ex- 
pressed by Eq. (2.5) is still valid in this case, provided 
we define the thermal potential in an appropriate way. 
We introduce 


h(x, y, 2, 0) = { cdé (8.1) 
0 
the total heat acquired by the unit volume. Further, 
we define a density function as 
F(x, y, 2, h) = = (8.2) 
0 0 


Finally, the thermal potential is defined for the volume 


v= ff fre 


Comparing with the definition in the linear case we see 
that F replaces the quantity (1/2)c#? which is obtained 
from Eq. (8.2) if c is independent of the temperature. 
With this definition we evaluate the variational 
quantity 


Tas 


(8.3) 


6V = SSS (OF /dh)éhdr (8.4) 
Conservation of energy requires 
h = —div H (8.5) 
Then, from Eq. (8.2), OF/Oh = 6 (8.6) 
Hence, 
= eaiv (6H )dr (8.7) 
which is formally identical with Eq. (2.6). The vari- 


ation 6D is defined in the same way as for the linear case. 
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We may then proceed exactly as in the linear case sy}. f 


stituting 61’ in Eq. (2.5) and deriving the differentiy 
equation of heat conduction from the variational prin, 
ciple. 

For the anisotropic nonlinear case the definition of | 
is the same as for the isotropic case and 6D remair; 
formally identical with the linear definition. The on] 
difference lies in the fact that the thermal resistivity 
matrix \;; is now a function of the temperature #. . 

The analogy of heat conduction to fluid flow in , 
porous material holds for the nonlinear case. Thy 
compressibility of the fluid and the permeability ar 
functions of the pressure. The mechanical potentia 
Il’ is then given by the same Eq. (8.3), where F is de 


F pdh (S.8 
0 


The incremental fluid pressure is p and / is the incre. 
mental mass of fluid acquired by the pores per unit 
volume of bulk material. It is related to the mass flow 
field H by relation (8.5). The definition of the dissi- 
pation function is identical with Eq. (3.10) and is pro- 
portional to the energy dissipated through friction. 

It should be noted that an equivalent formulation oj 
the variational principle is obtained in the isotropic 
case when the heat conductivity is a function only oi 
the temperature. By a well-known transformation we 
introduce the variable 


k(0)d6 (8.9) 
0 


fined as 


We may write 


k (00/Ox) = (Ou/Ox), etc. (8.10) 
c(06/Ot) = (c/k) (Ou/Ot) (8.11 

Eq. (2.1) then reads 
V7u = (c/k) (Ou/dt) (8.12) 


The problem is then mathematically identical with the 
case of a temperature distribution “ in a medium of 
thermal conductivity equal to unity and a heat capacity 
c/k function of u. 

The extension of the variational principle and the 
associated differential equations (3.14) for the general- 
ized coordinates is completely general. It includes 
the case where the surface heat-transfer coefficient A 
is a function of both the outside temperature and the 
wall temperature. Moreover, it may also be a given 
function of time. This feature leads to an exact 
treatment of problems which involve boundary-layer 
heat transfer under variable flow conditions and in- 
cluding radiation losses in the high temperature range. 
This is in addition to the inclusion of temperature de- 
pendence of heat capacity and thermal conductivity 
of the materials in cases where they vary widely in the 
considered temperature range. 

Surface radiation effects may be incorporated in the 
general expression for the dissipation function. We 
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may visualize a heat flow bifurcation at the surface. 
The heat flowing into the surface per unit area may be 
written /7, — H, when H, is the flow through the 
boundary layer and //, the heat flow lost by radiation. 
The radiation is equivalent to an additional branch of 
the thermal flow network. The radiation loss may be 
written 


H, = eo((7 + 0)? — T*] (8.13) 


where / + @ is the absolute temperature of the surface, 
¢ the emissivity, and o the Stefan constant. This is 


equivalent to 
H, = (8.14) 
with a temperature-dependent heat-transfer coefficient 
K, = (eo /6) [((7 + — 7" ] (8.15) 


If we wish to include radiation, the dissipation function 


(3.21) is replaced by 
D= dr + K, + 
l 
(8.16) 


where A, is the boundary-layer heat-transfer coeffi- 
cient. Note that the use of a reference temperature 
] is purely a matter of convenience. We could very 
well have put 7° = 0, in which case @ becomes the abso- 


lute temperature. 


(9) HEATING AND COOLING OF A SLAB WITH 
TEMPERATURE DEPENDENT PARAMETERS 


As an application of the method to a nonlinear 
problem let us consider a homogeneous slab with heat 
capacity and conductivity dependent on the temper- 
ature. We have shown that the problem may be re- 
duced to one of constant conductivity. Therefore, 
we shall consider the case where the heat capacity alone 
is a function of the temperature. This being a non- 
linear problem, different results will be obtained for 
heating and cooling. 

We first assume that the face at y = 0 is brought 
suddenly to a higher temperature 6) and that the heat 
capacity varies by a factor of two within the range of 
temperatures considered. 

= co[l + (0/0) ] (9.1) 
In Section (7) we have seen that there is a first phase 
where the temperature is approximately the same as 
in an infinitely thick slab. Using again the parabolic 
approximation we put 


6/0 = [1 — (v/a)? (9.2) 


where ¢, is the depth of penetration of the temperature 
rise. Following the general procedure of Section (8) 
lor the nonlinear case we evaluate 


[cao = + (Co/2) (07/00) (9.3) 
0 


F= coao = + (0/3) (67/0) (9.4) 
0 
The thermal potential is 
Vy = Fdy = (31/210) (9.5) 
0 


The heat flow // is obtained from 


H = (9.6) 


Putting ¢=1— (y/q) (9.7) 


we find 
H=q | = + (9.8) 
0 


Since Waa = — (@i/q) (9.9) 


the dissipation function is 


in 
D=— = af = 
2k Jo 2k 0 


0.0648... 
OK Co (9.10) 


The thermal force is obtained from 6// at y = O 
= = (13/30) (9.11) 
Hence, QO; = (13/30) (9.12) 
The differential equation for gq, is 
(31/210)coOo? + (0.0648, = 
(13/30)co%" (9.13) 


The solution with initial conditions g, = 0 at ¢ = Ois 


= 2.97 Viet (9.14) 


The transit time is found by putting q, = /, 
ty = 0.106 (9.15) 


We notice that the transit time for the present case 
where the heat capacity varies in the ratio of two to one 
is obtained by putting c = 1.2Sco in Eq. (7.14) for the 
linear case. This is appreciably lower than the average 
value c = 1.5c, indicating that the heat propagation 1s 
controlled more by the value of c in the region of lower 
temperature. 

It is of interest for comparison to consider the prob- 
lem of cooling instead of heating. In this case the face 
of the wall at y = 0 is cooled suddenly to a temperature 
6 = —@ at time ¢ = 0. We assume again a variation 
of the heat capacity from c = ¢) to c = 2c») between the 
extremes of temperature. This is expressed by the 
formula 


c = oll + (1/2) (6/40) ] (9.16) 


Because the problem is nonlinear the cooling problem 


= 


FLANGE 
A L 
~ 


Fic. 3. Cross section of wing structure. 


is different from that of heating even if it occurs between 
the same temperature limits. 
Proceeding exactly as above, we find 


V = (16/105) (9.17) 
H = /20)e — (1/3)¢*] (9.18) 
D = 0.0634 (co2/k)O0°qiGi2 (9.19) 
The thermal force is found by forming 6/7 at y = O or 
f=1; 
= = (17/30) (9.20) 
Hence the equation 
= 6.5(k/co) (9.21) 
or a, = 2.55V kt/co (9.22) 


The transit time ¢;—1.e., the time at which the opposite 
yall begins to cool off— is 


ty = 0.154(co/k)/? (9.23) 


If we compare this with Eq. (9.15) for the heating 
problem, we notice that it is quite a bit larger. In 
other words, the cooling takes about 40 per cent 
longer. Further comparison with the value (7.14) for 
the linear problem shows that the heat capacity be- 
haves as if it had a constant value 1.74 co which is higher 
than the average value 1.5 co. Therefore, the cooling 
occurs as if the heat capacity were higher than the 
average value while in the heating it is lower. In 
either case, the effective heat capacity tends to be 
closer to its value at the initial wall temperature.* 


(10) APPLICATION TO SUPERSONIC WING STRUCTURES 


We consider the supersonic wing structure whose 
cross section is shown in Fig. 3. It is heated on both 
sides by air brought suddenly at a temperature 4 
above the reference level 6 = 0 at time ¢ = 0. The 
heat is transmitted through the boundary layer to the 
flange and the web. We propose to calculate the tem- 
perature field in this structure as a function of time. 


* In general, we may conclude that the effective heat capacity 
will correspond to a temperature about midway between the 
average and the initial values. 
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In order to shorten the present calculation, certajy 
simplifying assumptions are introduced. The simplj 


fications are by no means essential to the methog 
However, they do lead to answers which are of accept. 
able accuracy in our problem. 

The present numerical solution must be considere 
as a first approximation which, although quite satis 
factory to our purpose, may be improved to includ 
any of the more complex features of the actual problens 
as pointed out in more detail below. 

We assume the heat flow to be one-dimensional—ig. 
the temperature is uniform across the thickness in both 
web and flange. This assumption is equivalent t 
stating that the transit time across the thickness is a smal] 
fraction of the total time history of the system. The 
transit time for a steel plate of 1/2 in. thickness is oj 
the order of 1 sec. The boundary-layer heat-transfer 
coefficient is taken to be a constant A. We consider 
the case where the wing is heated equally from both 
sides so that the problem is symmetric and no heat 
flow occurs at point .1/ middle of the web. The pro- 
cedure is easily adapted to the case of unsymmetric 
heating. The equivalent one-dimensional system js 
shown in Fig. 4. 

A first step is to calculate the temperature history for 
the flange alone in the absence of any web. If a de- 
notes the flange thickness and c its heat capacity per 
unit volume, the temperature of the flange #2 obeys the 
differential equation 

ac 6, = K(@) — 62) (10.1) 
The solution with the initial condition @. = 0 is 
6. = O[1 — (10.2) 
We have introduced the relaxation time of the flange 
boundary-layer system 
r =ac/K (10.3) 
In the numerical case considered below this relaxation 
time is found to be r = S6 sec. We now introduce the 
influence of the web. In the one-dimensional system 
of Fig. +, the web of actual thickness 2a, is represented 
by a plate of half the thickness, a;, attached end-on to 
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Fic. 4. Equivalent one-dimensional system with approxi- 
mate temperature distribution @ and heat inflow 1, through the 
boundary layer. 
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, plate representing the flange. In the first phase of 
the thermal history the heat has not yet had time to 
reach the mid-web 7. This phase occurs when the 
elapsed time ¢ is smaller than the transit time /, asso- 
ciated with the half length / of the web. This transit 
time according to Eq. (7.14) is 


t; = 0.0885(/?/«) (10.4) 


where x is the diffusivity of the web. In the numerical 
example below the half length of the steel web is 4.69 in. 
and the transit time of the half web is tf; = 105 sec. 

In the first phase we take advantage of a general 
feature exhibited by the approximate solution developed 
in Section (7) for the penetration of heat in a slab. 

We shall assume that in the web the penetration of 
heat is of a similar nature and that the temperature 
distribution in the web begins to rise at a moving point 
P, whose distance from the joint A of web and flange is 
determined by Eq. (7.12) for the penetration depth— 


1.€., 
P\A = q = 3.36V ct (10.5) 


Furthermore, we assume that the distribution of tem- 
perature in the web between the moving point ?; and 
the joint A is parabolic. This temperature in the web 


1S 
= 6,(y*/q") (10.6) 


where 0, is the temperature at point A and y is counted 
from P; as origin. Similarly, we assume that the in- 
fluence of the web on the temperature in the flange 
is to produce a drop below the temperature 62 calcu- 
lated above for the isolated flange, and that this drop 
penetrates to a point P whose distance from A is deter- 
mined by the same Eq. (10.5) for the penetration 
depth—i.e., PA = g. The temperature distribution 
in the flange between P and A is again assumed para- 
bolic. We write for this temperature 


0 = 6. + (0; — 62) (y?/q?) (10.7) 


where y is counted from P as origin. Finally, we also 
assume a parabolic distribution for the total heat /7/, 
per unit length which has flowed through the boundary 
layer between points A and P. We write 


H, = Hz + (Hs; — He) (y?/q*) (10.8) 


The quantity J/, represents the heat which has flowed 
into the flange without the web—i.e., outside of the 
region AP. Itis given by 


Fy = (10.9) 


It is a known function of time through Eq. (10.2) for 
#. The unknown quantities in the above expressions 
are H/;, which is the value of H, at the joint A, and 4. 
However, these two quantities are related by the law of 
conservation of total heat. We state that the total 
heat, which has flowed through the boundary layer be- 
tween A and P), is equal to the amount of heat stored 
in the web and flange between P; and P. This is ex- 
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pressed by 


A A A 
f H,dy = ac + ac f édy (10.10) 
P Pi 


Substituting the values (10.6) and (10.7) for @ and (10.8) 
for H, gives H; as a function of 4), 


Ai; = acBe, (10.1 1) 
with B=1+y, y=a/a 


We are thus left with a single unknown, the temperature 
6, of the joint. A differential equation for this quan- 
tity may be found by applying the general principles 
developed in the previous sections. The dissipation 
function of the boundary layer is 


. 
D = ( 2K) f H,*dy = (1/2K)H2*q + (1/3) X 
0 


(Hs — Hq + (1/10K)(Hs — (10.12) 
Substituting Eq. (10.9) and Eq. (10.11) yields 
D = (a*c?/K)q[(1/2)6.2 + (1 3)6.(86, — 6.) + 
(1/10) (86, — 6)?] (10.13) 


The generalized force 0; associated with the coordinate 
6, is found by evaluating the virtual work of the tem- 
perature on both sides of the boundary layer for a vir- 
tual change of heat flow. From Eqs. (10.8) and (10.11) 
we derive 


6H, = (v?/q?)6Hs = (10.14) 


The thermal force 0; is defined by 


A 
0:60; = (09 — (10.15) 


where 6 is given by Eq. (10.7). Hence, 
0, = (1/3)acBq(@ — 62) — (1/5)ac8q(@; — 62) (10.16) 
The differential equation for 4 is 
oD /d6, = (10.17) 
or 
(ac/3K)0. + (1/5K)ac(B6, — 6.) = 
(1/3) — 62) (1/5) (@; — @) (10.18) 


There are cancellations of terms due to Eq. (10.1). 
Moreover, by putting % -- 6; = 2 and introducing + 
from Eq. (10.3), Eq. (10.18) reduces to 


s+ Br = yrbp (10.19) 


With the initial condition z = 0 at ¢ = O and the value 
(10.2) for 42, the solution of this equation is 


z= — (10.20) 
This may be written 


6; = — (10.21) 


At this point it is of interest to introduce numerical 
values. The dimensions are taken to be 


| 
7 
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a ae 1 use the approximate value (10.21) of 6; to compute th, 
Oo ~ AIR TEMPERATURE time history of the temperature at any point in the we 
| by applying the Duhamel integral to the approxima, 
solution developed in Section (7). For instance, let 
Ag calculate the temperature at the mid-web .]/ in thy 
400 manner. According to Eqs. (7.25) and (7.26) a sudde 
unit rise of temperature at A at time ¢ = 0 produce 
no rise of temperature at the mid-web until ¢ >; 
After that at time = + Aft the temperature is 
(10.2 
200 + -- where 6 = 0.214/4. The transit time ¢, for the hal 
web as given by Eq. (7.14) is 
ty = 0.0885(I2/K) = 105 sec. 
By Duhamel's integral the temperature 6; at \/ is a 
time ¢; + Af 
0 100 200 300 400 500 600 Mt 
6; = [1 — e (10.23 
Fic. 5. Temperature 6; at the joint A and 6; at the midweb J/. 0 
Introducing the value (10.21) for 4, we find 
/ = 4.69 in., half length of web 
’ : D5 D5 — Al/Brj 9 
a = 0.375 in., thickness of flange ~ | 


2a, = 0.1 in., thickness of web 


The material is steel, for which 


K k/c = 0.0186 in.?/sec. 
¢ = 68:6 °F. 


diffusivity 
heat capacity per cu.ft. 


The heat-transfer coefficient of the boundary layer is 
assumed to be 


K = 90 B.t.u./hour ft.? °F. 


With these values the relaxation time of the flange- 
boundary-layer system is 


rT = ac/xk = 86 sec. 


The air is assumed to rise at ¢ = 0 to a temperature 
6) = 540°F. above the initial level. We introduce the 
value 8 = 1 + (a,/a) = 1.133. The numerical value 
of the joint temperature 6, as a function of time calcu- 
lated from (10.21) is plotted in Fig. 5. The present 
numerical values are the same as in the example treated 
by Pohle and Oliver (see reference 6),* who solved the 
heat conduction problem in the one-dimensional model 
by the exact method using the series and Laplace trans- 
form solutions of the corresponding partial differential 
equations. The procedure is quite elaborate. The 
exact value of 6; taken from reference 6 is also plotted 
for comparison in Fig. 5. It is seen to be in excellent 
agreement with our approximation. 

The result also points to the fact, which could have 
been surmised, that the influence on the temperature 
6, of the boundary condition at the mid-web J/ is 
not preponderant in this case. The assumption intro- 
duced for the first phase (¢ < f)) is practically valid 
for the complete time history of #;. We may therefore 


* The transient heat flow for the same geometry was also in- 
vestigated by Hoff,’ Parkes,* and Schuh.? 


This is plotted for times ¢ = ¢; + Afin Fig. 5 and con- 
pared with the exact value from reference 6. Again 
the two plots are in excellent agreement. 

Before terminating this example we shall briefly 
outline how the present method may be used to improve 
the accuracy of the solution and to include more com- 
plex features and refinements. 

A first step in this direction is to use the temperature 
history of the web calculated by expressions such as 
Eq. (10.23) instead of the parabolic distribution assumed 
above. The whole procedure is then repeated leading 
to an improved differential equation for 6;. A further 
improvement would take into account the calculated 
temperatures in the flange. This will lead, in general, 
to a differential equation with variable coefficients. 

Another improvement would be to get rid of the as- 
sumption of one-dimensional flow. For instance, we 
would start with the isolated flange problem and intro- 
duce instead of 6 two unknown temperatures, one on 
the side of the boundary layer, the other on the inside, 
with a parabolic temperature distribution across the 
thickness. We may then proceed further along similar 
lines. Two-dimensional flow in the web may also be 
introduced by a parabolic correction term for the tem- 
perature across the thickness along with a correction 
term for the nonuniformity of the heat flow near the 
joint. 

In the present example we have not mentioned the 
influence of the distance between webs. This effect 
may, of course, be taken into account by separating 
the problem into a third phase corresponding to the 
transit time associated with the distance between the 


webs. 

Antisymmetric heat flow results in the temperature 
A parabolic approximation may be 
It will tend exponen- 
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tially to the straight line distribution corresponding to 
steady-state heat flow along the web. 

Leakage from radiation and convection is also easily 
corrected for by introducing additional coordinates 
corresponding to the approximate leakage flow dis- 


tribution. 
Nonlinear problems, as we have seen by the example 


| of Section (9), may be handled by the present methods. 


For instance, if in the present example we wish to take 
into account the dependence of heat capacity and con- 
duction on temperature, we could first assume that 
points ? and P; move at speeds determined by approxi- 
mate nonlinear solutions obtained as in Section (9). 
Nonlinear leakage due to high temperature radiation 
may also be introduced. 

In many cases the heat-transfer coefficients of the 
boundary layer will be a function of both the time and 
temperature. This again imposes no restriction on the 
application of the method, the difference being simply 
that the differential equations in such a case will gen- 
erally have nonconstant coefficients. 

It may happen in such structures as we have con- 
sidered here that we must take into account the effect 
of a thermal resistance between the flange and the 
web. In such a case we should introduce the tem- 
peratures 6, and 6; of each side of the resistance. Pro- 
ceeding as above leads to one differential equation for 
these two unknowns. Another equation is obtained 
from the relation between the temperature drop through 
the resistance and the total rate of heat flow. We 
write 

— 63) = (1/3)e(d/dt) (64g) (10.25) 


where A, is the heat-transfer coefficient of the resistance. 
Finally we wish to call attention to the generality 
of the method of successive approximation illustrated 
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by the present example. We have just calculated the 
temperature 4, in the flange without the web and then 
introduced the influence of the web. In general it will 
be possible to proceed in a similar way by first calcu- 
lating the temperature fields in simple parts of the 
structure and introduce gradually various complexities 
and refinements. Each step will involve the solution 
of relatively simple first-order ordinary differential 
equations. 
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SUMMARY 


The optimum burning program for the horizontal flight of a 
rocket-powered aircraft is analyzed and discussed. The results 
of previous theories are extended to cover the general problem 
of minimizing an arbitrary function of the final values of time, 
mass, distance, and velocity. 

By using the indirect methods of the calculus of variations, it 
is shown that the totality of extremal arcs is composed of zero 
thrust sub-ares, sub-ares to be flown with maximum engine out- 
put, and variable thrust sub-ares. For the latter, a closed 
solution is obtained. 

Particular problems, such as maximum range, maximum en- 
durance, minimum propellant consumption, and maximum ve- 
locity increase, are treated within the general frame of the 
present theory for various types of boundary conditions. 

Detailed attention is devoted to the thrust programing which 
maximizes the range for the case where the overall flying time is 
prescribed. A method for constructing extremal paths is supplied 
under the assumption of a parabolic drag polar having either con- 
stant coefficients or coefficients depending on the Mach Number. 

An important difficulty associated with the linear aspect of the 
present problem is constituted by the fact that the Legendre- 
Clebsch condition fails to yield any information on the minimum 

or maximum character of the Eulerian paths. Moreover, the 
Weierstrassian function is zero at all points of the variable-thrust 
sub-are. This difficulty is overcome with a generalization of a 
previous method of the same author, based on the use of Green’s 


theorem. 
SYMBOLS 
a = speed of sound (ft. sec.~!) 
A = zero-lift drag (lbs. ) 
B = ratio of induced drag to square of mass 
(Ib. ~'ft.2sec. ~4) 
Cc = integration constant 
Coo = zero-lift drag coefficient 
D = drag (Ibs. ) 
dD, = partial derivative of drag with respect 
to velocity (lb. ft.~! sec.) 
pe = partial derivative of drag with respect to 
mass (ft.sec. ~?) 
E = maximum lift-drag ratio 
3 
F = 2, A\x/Jx = fundamental function 
k=1 
g = acceleration of gravity (ft.sec. ~?) 
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An Extension of the Theory of the Optimun 
Burning Program tor the Level Flight of a 
Rocket-Powered Aircraft’ 


ANGELO MIELE** 


Purdue University 


G; = function to be minimized 
H = constant defined by Eq. (36) 
A = particular value of the constant H 
Jx(K =1,2,3) = left-hand members of the equations rep: 
senting the constraints of the variation: 
problem 
K = ratio of induced drag coefficient to squx 
of lift coefficient 
= lift (Ibs. ) 
m = instantaneous mass (lb.ft.~'sec.?) 
my = reference mass (Ib.ft.~'!sec.?) 
Ms = reference mass (Ib.ft.~'sec.?) 
M = Mach Number 
p = atmospheric pressure (lb.ft.~?) 
= reference surface (ft.*) 
t = time (sec.) 
= thrust (Ibs.) 
u = V/V, = ratio of aircraft velocity to ex 
velocity 
V = absolute velocity of the aircraft (ft.sec 
V. = equivalent exit velocity of the rocket e 
gine (ft.sec.~!) 
X = horizontal distance (ft.) 
a = independent parameter of the rocket e& 
gine 
B = engine mass flow (lb.ft. ~'sec. ) 
“ = ratio of specific heat at constant pressur 
to specific heat at constant volume = 
1.4 for air 
A\k(A = 1,2,3) = Lagrange multipliers 
p = air density (lb.sec.*ft.~‘) 
® = function defined by Eq. (75) 
®, = function defined by Eq. (47) 
®, = function defined by Eq. (48) 
Vv = function defined by Eq. (76) 
VW; = function defined by Eq. (47) 
WV = function defined by Eq. (48) 
w = function defined by Eq. (63) 
1 = function defined by Eq. (61) 
we = function defined by Eq. (62) 
Subscripts 
i = initial condition 
f = final condition 
Superscripts 
* = derivative with respect to time 


= derivative with respect to Mach Number 


(1) INTRODUCTION 


ATTENTION has been devoted in re 
cent years to the theory of the optimization of the 
performance of aircraft and missiles. The classical 
quasi-steady approach, widely used in the performance 
analysis of conventional aircraft, has been frequently 
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repl 
take into consideration the inertia terms. 

When the acceleration terms are accounted for, sub- 
stantial modifications must be introduced into the 
analysis of the optimum paths. As a first approxima- 
ion. it may be said that—because of the inertia 
terms-—the nucleus of the problems of performance op- 
timization shifts from the domain of the ordinary theory 
of maxima and minima into the realm of the calculus 
of variations. This, in turn, implies a considerable in- 
crease in analytical difficulty, insofar as the optimum 
paths are no longer governed by equations in finite 
terms, but rather controlled by nonlinear systems of dif- 
ferential equations of higher order, integrable (in gen- 
eral) only by approximate methods. 

In connection with the applications of the calculus 
of variations to the general theory of aircraft and 
missiles performance, the fundamental problem is to 
determine both the geometry of the flight path and the 
mode of operation of the engine which optimize a given 
performance. Although for this problem general equa- 
tions have been indicated in the literature, neverthe- 
less the achievement of solutions in a closed form is 
practically impossible. It is for this reason that the 
attention of most researchers has been attracted by com- 
paratively simpler questions. In this connection, two 
main categories of variational problems have been in- 
vestigated. 

The first category (path programing problems) deals 
with the optimization of the geometry of the flight path 
for a prescribed mode of operation of the engine. 
Typical problems are the optimum climbing technique 
of a turbojet aircraft for a given power setting; or the 
optimum climbing technique of a rocket-powered air- 
craft for a given engine mass flow. 

The second category (thrust programing problems) 
deals with the optimization of the mode of operation of 
the engine for the case where the geometry of the flight 
path or an equivalent condition is given. A particu- 
lar case of this category of questions is the variational 
problem analyzed in the present article. * 


(1.1) Actual Status of the Research 


Previous investigations on the theory of the optimum 
burning program for a rocket-powered aircraft have 
been carried out by Hamel! and Tsien and Evans? 
for vertical flight and by Hibbs* and Cicala and Miele* 
for horizontal flight. The latter are reviewed below. 

The problem of maximizing the range in level flight 
for given end velocities and given propellant consump- 
tion, the flying time being free of choice, was first investi- 
gated by Hibbs.* He considered the ideal case of a 
rocket engine capable of delivering all mass flows be- 
tween zero and infinity and of an aircraft whose drag 
polar is parabolic with constant coefficients. Hibbs’ 
variational solution, however, did not satisfy the 
boundary conditions of the fixed end-points problem. 


*The present paper is a condensed form of the Technical 
Note 56-302 of the Air Force Office of Scientific Research.? 


As a consequence, he postulated the existence of two 
terminal ares to be traveled either in coasting flight or 
with a pulse-burning technique, and verified the above 
correct intuition by means of successive numerical 
analyses. 

More recently the same problem was attacked in 
reference 4 under more general hypotheses. Any re- 
striction concerning the shape of the drag polar was 
lifted. The more realistic case of an engine having a 
finite maximum burning rate was considered. For the 
case where a linear relationship between thrust and mass 
flow is assumed, a complete proof of the necessary and 
sufficient conditions for maximum range was obtained 
by Miele, using Green’s theorem. Concerning the 
more general nonlinear case, the proof was supplied by 
Cicala, by using the indirect methods of the calculus of 
variations and the new concept of index value. 


(1.2) Description of the Present Research 


Under the sponsorship of the Office of Scientific Re- 
search of the United States Air Force, an investigation 
of a wide class of variational problems of interest in the 
theory of the burning program has been undertaken. 
The final objective is an attempt to generalize all the 
existing thrust programing theories. 

Within the general frame of the aforementioned re- 
search, the burning program for level flight has been 
reexamined with the following general problem in 
mind. “‘A rocket engine, capable of delivering a vari- 
able thrust, is applied to a vehicle of given configura- 
tion moving along a straight, horizontal path. The 
initial time instant (f;), mass (m;), abscissa (Y,), and 
velocity (1’;) are prescribed; some, but not all, of the 
final coordinates m,, X,, are given. It is required 
to determine the thrust-programing technique which 
minimizes a specified, though arbitrary, function 
Gy (ty, my, Xs, Vy) of the coordinates of the final point.”’ 

Particular cases of the above general problem are 
the following: (1) the case G; = —Xy,, maximizing 
the range; (2) the case Gy; = —/,, maximizing the time; 
(3) the case G; = —my, minimizing the propellant 
consumption; and (4) the case G, = —I’,, maximizing 
the final velocity. 

The most interesting by-product of the above anal- 
ysis is, perhaps, the investigation of the burning pro- 
gram which maximizes the range X, for the case where 
the initial and final values of velocity, mass, and time 
are prescribed. The practical importance of this prob- 
lem was indicated at the September 1955 meeting of 
the American Rocket Society by Mr. Jack Lorell of the 
Jet Propulsion Laboratory, California Institute of 
Technology. In the discussion following the presenta- 
tion of the paper of reference 4, Mr. Lorell inquired 
about the possibility of extending the analytical tech- 
niques developed by Professor Cicala and the writer 
to the case where a maximum range is desired for a 
prescribed value of the time. He pointed out that 
paths of maximum absolute range (free time, see refer- 
ences 3 and 4) are relatively slow, so that it is of engi- 
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Parametric representation of the engine characteristics. 


neering interest to study the burning program for the 
case where a specified value is prescribed for the overall 
time of flight. 


(2) THE Puysics OF THE PROBLEM 


In the following sections of this paper, these hy- 
(a) the rocket-powered aircraft is 
ideally regarded as a particle of mass (m) variable with 
the time (f); (b) a rectilinear level path is considered; 
(c) the small angle between the vectors thrust (JT) and 
velocity (V) is neglected; (d) the equivalent exit 
velocity (V’,) of the rocket engine is regarded as a con- 
stant, independent of the operational condition of the 
engine; (e) the engine is capable of delivering all mass 
flows (8) bounded between a lower value (6 = 0) and an 
upper value (8 = Bnar); and (f) the aerodynamic lag 
is disregarded—i.e., lift (L) and drag (D) forces are 
calculated as in unaccelerated flight. 

Because of the hypotheses (a), (b), and (c), the 
equations of motion on the tangent and on the normal 
to the flight path are written as follows: 


potheses are used: 


T—-D-—-mV=0 
L—mg=0 


where g is the acceleration of gravity and the dot sign 
denotes derivative with respect to time. 


(2.1) The Drag Function and Its Derivatives 


In view of the assumptions (b) and (f) the drag D 
can be considered as a function D (V’, L) of velocity and 
lift, even accounting for the effects of the Mach and 
Reynolds Numbers on the polar. 
L can be eliminated by means of Eq. (2) and the drag 
can be reduced to a function of velocity and mass only. 
In the following discussion, two cases are considered: 
(a) the general case of a drag function of the form 


In turn, the lift 


D = D(V, m) 
and (b) the particular case of a parabolic drag polar 


D= A+ Bm? 
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where 
A (1 2) CpopS V? (5) 
B = 2Kg?/pSV? (6 
In the above equations Cp» is the zero-lift drag coef}. 
cient, A the ratio of the so-called induced drag coef}. 
cient to square of lift coefficient, S a reference surfag 
and p the air density. From Eq. (4) the partial deriva. 
tive of drag with respect to mass can be calculated as 


[OD, Ont |ymconst. = 2Bm (7 


Then Eqs. (4) to (6) yield the partial derivative of drag 


with respect to velocity, 


[2 + + (Bm?/V) [—2 + M(K'/K)| 


(8) 

where 
Ca’ = dC po dM (9) 
= dK/dM (10) 


It is to be noted that Eq. (8) has been calculated by 
neglecting the variation of drag associated with vari- 
ations of the Reynolds Number—-i.e., by considering 
Cy» and K as functions of the Mach Number (.1/) only. 


(2.2) The Thrust Function 


Because of hypothesis (d) the thrust 7 is expressed as 
a linear function of the mass flow @: 


T = BV, (11) 
In turn, the hypothesis (e) implies that the inequality 
BS Braz (12) 


be satisfied by the class of all arcs investigated in the 
present report. 

The limitation (12) which the size of the engine im- 
poses on the mass flow 8 has a critical importance in 
determining the behavior of the solution, as has been 
recognized in references 3 and 4. In reference 4 the 
above inequality has been handled by using the concept 
of index value. Another possibility—-perhaps more 
elegant from a formal point of view—consists of using a 
parametric representation of the engine characteristics. 
Such a technique is developed in the present report. 

The mass flow @ is represented as a function of a 
parameter a having the following properties: (a) 
for —*” < a@< a, the mass flow is 8 = 0; (b) for a < 
a < + o, the mass flow is 8 = Bngz; and (c) for a < 
a < a, a one-to-one correspondence is assumed be- 
tween 6 and a. 

With the above scheme, a is considered as the inde- 
pendent parameter of the rocket engine and is allowed to 
vary between —o© and +. The mass flow £ be- 
comes a dependent quantity, varying between 0 and 
Bmax according to the scheme of Fig. 1. In turn, the 
thrust 7 is limited between 0 and 7,,,., according to 
Eq. (11) and Fig. 1. 

Notice that the condition d8/da = 0 represents 
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either a coasting flight or a flight with maximum engine 
output. On the other hand, d8/da # 0 represents any 
other operating condition of the engine intermediate 
between the two limiting ones. Notice also that a 
is only a parameter and that there is no necessity of 
attributing to it any special physical meaning. 


(3) VARIATIONAL FORMULATION 
The following set of equations* is now considered: 
J; =m + B(a) = 0 (13) 
Jo =X-V=O (14) 
Js = V+ (DV, m) — VeB(a)]/m} = 0 (15) 


Eq. (13) is the relationship between time-rate of 
change of the airplane mass and engine mass flow 8. 
In turn, 8 is a known function of a, according to the 
parametric representation illustrated in Fig. 1. Eq. 
(14) is the differential relationship between space and 
time. Eq. (15) is derived from the equation of motion 
on the tangent to the flight path after accounting for 
the analytical expressions for thrust and drag. 

The time (¢) is assumed as the independent variable 
of the problem. The dependent variables are the mass 
(m), the space flown horizontally (X), the velocity 
(V’), and the engine parameter (a). 

After prescribing the initial coordinates f; = 0, m; = 
m(0), X; = 0, V; = V(O) and some, but not all, of the 
coordinates at the final point, the variational problem 
is stated as follows. ‘Between all sets of functions 
m(t), X(t), V(t), a(t) satisfying Eqs. (13) to (15) and 
the prescribed end-conditions to determinej the special 
set which minimizes a specified function G,(t,, my,, X,, 
l’,) of the coordinates of the final point.”’ 

It is to be noted that no condition can be imposed on 
a, and ay, i.e., on the mass flows 8, and 8, In view 
of the analytical nature of the problem, the end-values 
for 6 are a mathematical consequence of the set of 
Euler equations, as shown in the following sections. 


(3.1) Euler Equations 


A set of variable Lagrange multiplers \,(/), A2(/), 
\;3(t) is introduced and the following expression formed: 


(16) 
1 


where J;, Jo, J3 denote, respectively, the first members 
of Eqs. (13) to (15). Since the unknown functions are 
four in number, one must write four Euler equations. 
These are written as 


(d/dt) = OF/dzy (J = 1, 2,3, 4) (17) 


* In the terminology of the calculus of variations, Eqs. (13) to 
(15) are nonholonomic or differential conditions. 

} In this way, the variational problem is formulated within 
the general frame of the problems of Mayer type. The Mayer 
problem is known as the most general problem of the one- 
dimensional calculus of variations,» equivalent in scope and range 
to the formulations of Lagrange and Bolza. 
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where 2; = m, 2 = X, 23 = V and 3% = a. From 
Eqs. (13) to (17) the following explicit form of the 
Euler equations is obtained: 


= — D)/ + (D,, m)} Xs (18) 


he = 0 (19) 
As = + As(Dy m) (20) 
= (dg da) [Ay V. m) | (21) 


(3.2) First Integral 


Since the time does not explicitly appear in the prob- 
lem, the first integral shown below is valid > 


4 
(OF /03,)3, = C (22) 
=1 


where C is an integration constant. After accounting 
for Eqs. (13) to (16), Eq. (22) is rewritten as follows: 


Aint + AX + ASV = C (23) 


The above expression is a mathematical consequence of 
Eqs. (18) to (21) and may be used to replace one of 
them-——for instance, Eq. (18). It must be emphasized 
that the set of Euler-Lagrange equations (18) to (21) 
and the first integral (22) or (23) hold for all problems 
investigated in this report, regardless of the nature of the 
function my, Xs, and of the end-conditions. 


(3.3) Discontinuity of the Eulerian Solution 

The Euler equation (21) is particularly important, 
because it shows that the extremal arc is discontinuous, 
being generally composed of (a) sub-ares of equation 


dB/da = 0 (24) 
and (b) sub-arcs of equation 
Ar — A3(V./m) = O (25) 


According to the parametric representation shown in 
Fig. 1, Eq. (24) represents either a coasting flight or a 
flight with maximum engine output. Eq. (25), on the 
other hand, represents a flight condition with a con- 
tinuously variable thrust, as shown below. 


(3.4) Erdmann-Weierstrass Corner Conditions 


In view of the discontinuous character of the solu- 
tion, the Erdmann-Weierstrass corner conditions must 
be applied. These are continuity conditions to be 
satisfied at every corner of the discontinuous extremal 
solution by the following five quantities: 

J=1 03, 

Since F is independent of a, the continuity of OF 0a 
is inherently verified at all corner points. The analo- 
gous requirement for the remaining four quantities 
(26) leads to 


(Ax)— = (Ax)+ (K = 1, 2, 3) (27) 
(C)_ = (C)+. (28) 
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where the subscript (—) denotes a condition immedi- 
ately before the junction and the subscript (+) a condi- 
tion immediately after the junction. The Erdmann- 
Weierstrass conditions require, therefore, (a) that the 
multipliers \x(A = 1, 2, 3) be continuous at the junc- 
tion points, and (b) that the value of the integration 
constant C be the same for all sub-ares of the discon- 
tinuous extremal solution [this latter condition implies, 
also, that the quantity \; — A3(V,/m) be zero at a junc- 
tion point; thus a corner point may lie on the variable- 
thrust sub-are defined by Eq. (25) ]. 


(4) SOLUTION OF THE PROBLEM IN THE MAsS-VELOCITY 
PLANE 


In the following section general solutions are pre- 
sented for predicting the optimum burning program in 
the mass-velocity plane. These solutions involve, in 
general, a number of constants to be determined accord- 
ing to the nature of the minimal problem and of the 
end-conditions. 


(4.1) Constant Thrust Sub-Arcs 


Because of Eq. (13), the acceleration |’ is rewritten 


as 
= —B(dV/dm) (29) 

so that Eq. (15) yields 
dm/m = [8/(D — BV,)]dV (30) 


Eq. (30) is a first-order differential equation involving 
V and m and must be integrated—in general—by ap- 
proximate methods, even if the drag polar has the para- 
bolic form represented by Eq. (4). Exact solutions 
are possible only in special cases, such as for instance 
the coasting flight (zero mass flow) where 


B= 0; 


and the pulse burning (infinite mass flow) where 


m = const. (31) 


8= mexp(V/V.) = const. (32) 


(4.2) Variable Thrust Sub-Arcs 


For the sub-ares of the extremal solution governed 
by Eq. (25) the Lagrange multipliers \; and \; can be 
eliminated with the following procedure. 

After combining Eqs. (13), (14), (15), (23), and (25), 
extensive manipulations yield 

(V. D) [Ao V 
As = (m/D) [eV — (34) 
By calculating the time-derivative of Eq. (34) and 
accounting for Eqs. (13), (15), (19), and (20), the fol- 
lowing result is obtained: 


((V/V.) H] (D + V.Dy — mD,) — D =0 
H = 


(33) 


(35) 
where (36) 
Eq. (35) is a relationship in finite terms involving mass 
and velocity which holds whatever be the shape of the 
function G; to be minimized. The constant 7 must be 
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determined, for each particular protlem, according { 
G,and the end-conditions. 

For the special case of a parabolic drag polar wit) 
coefficients depending on the Mach Number, Eqs. (4 
to (10) and (35) yield an explicit solution for the mas. 
velocity relationship 


m/m, = M? V (Fy — HF.)/(F; — HF.) (3 
where 
m, = ypS/ 2g (38 
Fy = Cpoo(1 + + WCppo’ (39 
F, = + 2/u) + Coo’ V,/a (40 
F; = K(3 + u) — MK’ (41 
FP, = K(1 + 2/u) — K’'V./a (42 


In the above equations m, is a reference mass, y the 
ratio of specific heat at constant pressure to specific 
heat at constant volume (air), a@ the speed of sound 
u = V/V,, and p the atmospheric pressure. For the 
sub-case where the coefficients of the parabolic polar 
are constant, Eq. (37) simplifies into 


m 8 
= 3) 
M» u+3 — A(1 + 2/n) 
where my» is another reference mass defined by 


mz = (pSV,2/2g) V Cmo/K (44) 


(4.3) Determination of Distances and Times 


The determination of distances and times in the (J 
m) plane is accomplished with the following formulas, 
deduced from Eqs. (13) to (15): 


f 
= [did V + VW, dm} (45) 
= [od + Wodm | (46) 
where 
®,(m, 1) = —(mV/D), V) = 
—(V.V/D) 
V) = —(m/D), Wom, V) = D) (48) 


(5) BouNDARY CONDITIONS 


For the problem under consideration, the boundary 
conditions include a number of fixed end-point condi- 
tions plus a number of natural conditions. The latter 
must be deduced from the following general trans- 
versality condition 


OF 
6G, + | dé, (6z7 — 28) | =) (49) 
which is to be identically satisfied for all systems of 
variations consistent with the prescribed end-condi- 
tions. After accounting for Eqs. (13) to (16), for the 
first integral (22) and for the fact that the initial co- 
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ordinates are given, Eq. (49) leads to 
6G; + + + V7 Cét ly = (50) 
Since G, depends on my, Xy, Vy, ty, Eq. (50) splits into 
the following four sub-conditions: 
Om,) + Ary] bm, = O (51) 
[(OG, OX do] 6X, = 0 
(OG, Ol’,) + Asy] = 0 


(6) SpecrAL PROBLEMS 


With the help of the generalized solutions discussed 
in Section (4) and of the boundary conditions indicated 
in Section (5), a number of particular problems are now 


analyzed. 
(6.1) Maximum Range, Free Time, Given Final Velocity, 
Given Propellant Mass* 


With reference to the problem of maximizing the 


range (G, = —X,) for a given final velocity (6V, = 0), 
given final mass (6m, = 0) and flying time free of choice 
(ot; = arbitrary), the transversality condition (51) 


vields: C = 0,4» = 1. For adrag function of the form 
D=D (I, m) the optimum burning program includes 
constant thrust sub-ares, defined by Eq. (30), and vari- 
able thrust sub-ares whose expression can be derived 
from Eq. (55), 


+ V.Dy — mD,,) — VD = 0 (52) 


For a parabolic polar with coefficients Cpo, A depending 
on the Mach Number only, Eq. (52) yields 


me. + u) + 


53 
K(3 + u) — MK’ 


my 
an equation which can also be determined as a particu- 
lar case of Eq. (37) for 7 = 0. In the case where the 
coefficients of the parabolic polar are constant, Eq. 


(53) simplifies into* 
= (1 + u)/(3 + n) (54) 


No further attention is given to this problem. As a 
matter of fact, the way in which the constant thrust 
sub-ares (8 = 0 or B = Bmar) must be combined with 
the variable thrust sub-ares (53) or (54) to satisfy the 
prescribed boundary conditions has been already thor- 
oughly discussed by Hibbs* and by Cicala and Miele.‘ 


(6.2) Maximum Endurance, Free Horizontal Distance, 

Given Final Velocity, Given Propellant Mass 

The problem of maximizing the time (G, = —t,) for 
a given final velocity (6, = 0), given final mass (6m, = 
0) and horizontal distance free of choice (6X, = arbi- 
trary) is now discussed. Though this problem has a 
minor technical importance, nevertheless it is of interest 
to remark that its solution is embodied in the general 


*The problem of maximizing the range for a given time is 
treated in detail in Section (9) of the present article. 
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theory presented here. The transversality condition 
(51) yields C = —1, \2 = O, so that Eq. (35) simplifies 
into 


D+ V.Dy — mD,, = 0 (55) 


For a parabolic drag polar with coefficients depending 
on the Mach Number only, Eq. (37) yields 


m |Cool(2 + u) + 
mM, K(2 + u) — MK 
which simplifies into 
= (57) 


if Cy) and K are constant. The solution is therefore 
composed of sub-ares satisfying Eq. (56) or Eq. (57), 
coasting flight (8 = 0), and sub-ares to be flown with 
maximum engine output (8 = Bye). To satisfy the 
boundary conditions of the problem, the above sub- 
arcs must be combined in analogy with the general pro- 
cedure discussed in reference 4 for the problem of Section 
(6.1). 

It must be emphasized that, for the two problems 
discussed in Sections (6.1) and (6.2), the optimum 
burning program can be predicted in the (1’, m) plane 
independently of distances and times. Once the opti- 
mum mass-velocity relationship is known, the values 
for XY, and ¢t, are found by subsequent evaluation of the 
integrals (45) and (46). 


(6.3) Minimum Propellant Consumption, Free Time, 

Given Range, Given Final Velocity 

The problem of minimizing the propellant consumed 
(G; = —m,) for a given final velocity (61, = 0), given 
range (6X, = 0), the fiying time being free of choice 
(6t; = arbitrary), is now analyzed. The transversality 
condition (51) yields \y, = 1, C = Oso that the expres- 
sion for the variable thrust sub-arcs reduces to Eq. 
(52). It is therefore concluded that the problem of 
maximizing the range for given propellant consumption 
and the problem of minimizing the propellant consump- 
tion for given range /ead to the same totality of extremals. 
Such a property is due to the fact that the system of 
Euler equations (18) to (21) is homogeneous in the 
Lagrange multipliers. As a consequence, two propor- 
tional sets of multipliers may be associated with the 
same solutiont m(t), X(t), V(t), a(t). 


(6.4) Maximum Final Velocity, Free Time, Given Range, 

Given Propellant Consumption 

For the problem of maximizing the final velocity 
(G; = —Vy), for a given distance (6X; = 0), given 
mass of propellant (6m; = 0) and free time (dt, = 
arbitrary) the transversality condition (51) yields 
\sy = 1,C = 0. This problem leads, therefore, to the 
same totality of extremals associated with the problems of 
Sections (6.1) and (6.3). 

7 The above property can also be regarded as a by-product of 
the general Mayer’s law of reciprocity of isoperimetric prob- 


lems.® 
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(6.5) Maximum Final Velocity, Free Time, Free Range, 
Given Propellant Consumption 

For a given propellant mass (6m, = 0) and for the 
problem G, = —V,, 6t; = arbitrary, 6X, = arbitrary, 
the condition (51) yields \3, = 1, C = 0,42 = 0. The 
variable thrust sub-are (35) is excluded from the solu- 
tion, which only consists of constant thrust trajectories. 
Obvious physical considerations rule out the branch 
8 = 0. As a consequence, the extremal are for the 
above problem degenerates into a single flight path 
to be flown with maximum engine output. 


(7) LEGENDRE-CLEBSCH AND WEIERSTRASS CONDITIONS 


According to reference 5, the necessary condition 
due to Legendre-Clebsch—in its strengthened form— 
requires that the following inequality be satisfied at 
all points of the Eulerian arc for G, to be a minimum: 
4 4 O°-F 
— > 0 (58) 

where 7, = ™, m2 = xX.n - V and m = a. The 
system of variations 6m, 6X, 61’, 6a must be consistent 
with the equations of motion—.e., with 


(K=1,2,3) (59) 


The development of Eq. (58) yields 


Q = (0°F/2a°) (6a)? = 
— As(V-/m)] (d*B/da*) (da)? (60) 


For the constant thrust sub-arcs, defined by Eq. (24), 
the second derivative d°8/da? is zero, so that Q = 0; in 
addition, for the variable thrust sub-are Q = 0, because 
of Eq. (25). Therefore, it is concluded that the Legen- 
dre-Clebsch condition does not yield any information 
on the minimum or maximum character of the Eulerian 
path. 

An analogous difficulty is encountered when applying 
the more stringent Weierstrass condition. As a matter 
of fact, the Weierstrassian function? is zero at all points 
of the variable-thrust sub-arc. 

The weaknesses indicated above have their origin in 
the fact that the present variational problem is linear 
in m, V, X, and 8. These difficulties, however, can 
be offset with the help of the properties of the three 
functions w, w:, and w defined in Section (8) and of 
Green's theorem as used in Section (9) of this report. 


(8) DEFINITION OF THE FUNCTIONS 1, w, AND w 


The following functions of the variables V and m 
are now introduced: 


wi(V, m) = (0W,/0V) — (04,/dm) (61) 
w(V, m) = (OW2/dV) — (O42/dm) (62) 
w(V, m, H) = — HVan (63) 


where V2 are defined by Eqs. (47) and (48) 
and by Eq. (36). 
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After accounting for Eqs. (47) and (48), Eqs. (61) t 
(63) yield 


wD? = V(D + V.Dy mD») V.D (64) 
= D + V.Dy mD,, (65 
wD? = (D + V.Dy mD,,) (V HV.) VD (66 


For a parabolic polar having coefficients C pp and K de. 
pending on the Mach Number only, Eqs. (64) to (66) 
in combination with Eqs. (4) to (10), yield 


wD? Cyo/AV, = Fy — (F3/M*) [m/m (67) 
weD?C po A = (F; M"*) [m m,|* (68) 


wD?Cpy/AV, = (Fi — HF.) + 
— [m/m,]? (69) 


where Fi, F2, F3, Fy are defined by Eqs. (39) to (42), 
If the coefficients of the parabolic polar are constant, 
Eqs. (67) to (69) lead to 


wD?/AV, = 1+u — [(3 + u)/u‘] (70) 
wD?/A = [1 + (2/u)] [1 — (1/u*) (m/m.)?] (71) 


4 
AV. u 
+ (2/u)] k | 
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A comparison between Eqs. (35) and (66), (37) and 
(69), and (43) and (72) shows that the points of the 
(lV, m) plane where w = 0 belong to the variable thrust 
sub-ares associated with all problems of the form G,(t, 
my, Xs, Vy). Egs. (35), (37), and (43), therefore, can 
also be designated as w(m, I’, H7) = 0. 


(9) ANALYSIS OF THE MAXIMUM RANGE PROBLEM FOR 
A GIVEN TIME 


For the problem of maximizing the range (G, = 
—X,) for given final velocity (6V, = 0), given final 
mass (6m, = 0) and given flying time (6¢; = 0) the trans- 
versality condition (51) yields \, = 1. 

The extremal arc is, therefore, composed of sub-ares 
to be flown with engine shut-off, sub-arcs to be flown 
with maximum engine output and variable thrust sub- 
ares, the latter being defined by Eq. (35) for a general 
drag function of the form D = D(V, m). The next 
step is to determine how the constant thrust sub-ares 
must be combined with the variable thrust sub-ares 
to satisfy the prescribed boundary conditions and to 


‘establish whether the Eulerian solution yields a mini- 


mum or a maximum for the range. 

The above topics can be substantially clarified if the 
discussion is divided into two cases; (I) the case where 
the coefficients Cyo, K of the parabolic polar are con- 
stant, and (II) the case where Cp) and K are functions 
of the Mach Number. Only the first of the above 


cases is considered in the present article for the sake of 
brevity. Concerning the second, the reader is referred 
to the extensive analyses carried out by the writer in 
reference 7. 
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9.1. The Case of a Parabolic Polar With Constant 

Coefficients 

For this highly idealized case the extremal arc is 
composed of constant thrust sub-ares (8 = 0 or 8 = 
8nar) and sub-ares satisfying Eq. (43). The latter can 
also be indicated as w{m, I’, /7) = O according to the 
discussion of Section (S). 

(9.1.1) Family of Variable-Thrust Sub-Arcs-—Eq. 
(43), which represents the family of variable-thrust 
sub-ares, is plotted in Fig. 2 for different values of the 
constant /7. 

For J] @ 0 Eq. (43) is such that the velocity is a 
single-valued function of the mass m. On the other 
hand, for /7 > 0 the velocity |” is a double -valued func- 
tion of the mass m, for such a case, therefore, each vari- 
able thrust sub-are splits into two branches—a _ high- 
speed branch (solid line of Fig. 2) and a low-speed 
branch (dotted line). 

(9.1.2) Minimizing or Maximizing Properties of the 
Variable-Thrust Sub-Arcs—-As it was pointed out in 
Section (7) both the Legendre-Clebsch condition and the 
Weierstrass condition fail to yield any information on 
the minimum or maximum nature of the Eulerian arc. 
In spite of this, important information can be obtained 
by using Green’s theorem in combination with the 
special properties of the functions @), w2, and w defined 
in Section (8). In this connection, the following re- 
marks are important. 

(a) The sub-arc J7 = O is the geometrical locus of 
the points where w (I, m) = 0. The function @(1’, 
m) is positive in the region A (Fig. 3) located on the 
right of wa, = 0 and negative in the region B located on 
the left of w; = 0. 


AX, — HV At, = [b(m, V, + v(m, V, — 


1.€., 
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(b) The sub-are J] = + © is constituted by the suc- 
cession of points where wo(I", m) = O. The function 
wo( 1", m) is positive in the region A (Fig. 3) located on 
the right of w. = 0 and negative in the region B located 
on the left of w. = 0. 

(c) For finite values of /7 such that /7 < 0 the funce- 
tion w(m, |’, H) is positive in the region A (Fig. 4) 
located on the right of w = 0 and negative in the region 
B located on the left of w = 0. 

(d) For finite positive values of //, the function 
w (m, I’, I) is positive in the region A (Fig. 4) located 
on the right of the solid line w = 0 and in the region C 
on the left of the dotted line w = 0. On the other 
hand, w is negative in the intermediate region B located 
between the solid line and the dotted line of Fig. 4. 

To establish whether the trajectories of Fig. 2 yield 
maximum or minimum solutions, a particular value for 
H is considered (for instance // < 0, Fig. 4) in connec- 
tion with the case where both end-points / and F are 
located on the variable thrust sub-are defined by the 
selected value for /7. The Eulerian sub-are /QF is 
compared with the control path* /PQRF and the fol- 
lowing differences are formed: 


AX, = [Xr lige LX (73) 
= — Itrliperr (74) 


After introducing the multiplying factor //1°, and the 
definitions 


V, = — (75) 
wv, — (76) 


Wim, 


one obtains 


V, + v(m, V, (77) 
IPORF 


AX, — HV At, = [i(m, V, + v(m, V, + [b(m, V, + v(m, V, (7s) 
IOPI OFROQ 


By the use of Green's theorem the two line integrals appearing in Eq. (7S) are transformed into two surface in- 
tegrals, respectively associated with the areas S, and S, internal to the contours /PQ/ and QFRQ: 


Ar ay = J s, Om I. 


ol 


w(m, V, W)dVdm + ff wim, V, H)dVdm (79) 


where the function w is defined by Eq. (63). The 
minus sign preceding the first surface integral is due to 
the clockwise sense of the circuit /?Q/. Analogously, 
the plus sign preceding the second surface integral 
stems from the counterclockwise character of the cir- 
cuit QVFRQ. Since the function w (m, I’, /7) is—for HT 


< 0 negative in the region S, and positive in the re- 
gion S,, the following result is obtained: 


OP 
a ) dVdm = 


om 


AX, — HV_.At,>0 (SO) 


In the case where the comparison path /PQRF and 
the Eulerian sub-arc /QF are flown in the same time 
(At, = 0), the inequality AY, > O follows. This means 
that the Eulerian are /QF maximizes the range for a 


* The comparison path 7PQRF must be consistent with the 
limitation 0 < 8 < Baz. 
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Fic. 2. 
0 for the case of a parabolic polar with constant coefficients. 
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Family of variable thrust sub-ares w(m, V, H) = 


‘w(mV,H) = 0 
m 
‘8 
w(m,V,H) <0 w(m,V,H) >O 
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H<O 
w=0 w=0 
/ 
© / ® 
wim,V,H) >0 wim,V,H) <0 w(m,V,H) 
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Fic. 4. Properties of the function w (m, V, H). 
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Fic. 6. Family of extremal ares corresponding to given initial 
and final values of mass and velocity (parabolic polar with con- 
stant coefficients). 
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Fic. 5. Example of combinations of sub-arcs. 
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prescribed flying time. On the other hand, if the com- 
parison path is chosen in such a way that AX, = O, the 
inequality At; > 0 follows, since 7 has been assumed 
negative. This means that the Eulerian path /QF of 
Fig. 4 maximizes the flying time for a prescribed range. 

By similar arguments it can, furthermore, be shown 
that: (a) the sub-are J/7 = 0 is associated with maxi- 
mum range solutions for the case where the flying time 
is free of choice (Fig. 3), (b) the sub-are T = + © is 
associated with maximum time solutions for the case 
where the range is free of choice (Fig. 3), (c) the high- 
speed branch of all sub-ares where /7 > 0 (solid line of 
Fig. 4) maximizes the range for a prescribed time and 
conversely—minimizes the time for a prescribed range, 
and (d) the low-speed branch of all sub-ares where // > 
() (dotted line of Fig. 4) minimizes the range for a pre- 
scribed time and—-conversely——maximizes the time for a 
prescribed range. 

(9.1.3) Simple Criterion for Recognizing the Character 
of the Variable-Thrust Sub-Arc w(m, V, H) = O—With 
reference to problems where the flying time is pre- 
scribed, the following criterion is derived from the 
discussion of Paragraph (9.1.2). “‘A variable thrust 
sub-are w (m, I’, HZ) = Ois associated with a maximum 
range solution if the function w is negative in a neigh- 
borhood immediately on the left of a = 0 and positive 
in a neighborhood immediately on the right; the same 
sub-are yields a minimum range solution in the case 
where the signs of the function w are reversed.”’ 

(9.1.4) Combination of Sub-Arcs—Consider one par- 
ticular sub-are of the family defined by Eq. (43)—for 
instance, the sub-arc /] = —0.2. Assume now that 
the initial point J and the final point F do not belong 
to the curve w (m, I’, H7) = 0. 

The problem is to determine how the above variable 
thrust sub-are w = 0 and the constant thrust trajectories 
8 = 8B = Byar Must be combined. It is clear that, 
in order to satisfy the boundary conditions associated 
with the prescribed end-values for mass and velocity, 
the Eulerian are connecting the points / and F must be 
generally composed of three branches (see Fig. 5): an 
initial branch (/.1/); a central branch (.1/.V) to be flown 
along w(m, V, 7) = 0; and a final branch (VF). 

The initial sub-are /.1/ is to be flown at the maximum 
burning rate of the engine if point / is located at the 
left of w = 0, or in coasting flight if point J is placed 
at the right of w = 0. The final sub-are NF is to be 
flown with engine shut-off if point F lies at the left of 
# = 0, or at the maximum burning rate of the rocket 
engine if F is located at the right of w = 0. In other 
words, four possible types of boundary conditions may 
exist‘ in the mass-velocity plane. They are designated 
with the Roman numerals I, II, III, and IV and are 
shown in Fig. 5. 

(9.1.5) Solution of the Boundary Value Problem 
Assume now that the seven quantities ¢;, m;, V;, 
ty, my, Vy are prescribed and that it is desired to maxi- 
mize the range X’y. 

Consider (Fig. 6) the family of paths /ABF, ICDF, 
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Fic. 7. Distance flown and overall time as a function of the 
parameter H (parabolic polar with constant coefficients). 


IEGF, IRJF, IKLF composed of sub-ares 8 = con- 
stant and sub-arcs w(m, V, 7) = 0. Each of the above 
paths is a solution of the set of Euler equations (1S) 
to (21). Moreover, all of them satisfy the prescribed 
initial and final conditions concerning mass and ve- 
locity and the Erdmann-Weierstrass corner conditions.’ 

The subsequent step is to account for the time con- 
straint. No direct solution seems possible, and the 
following indirect procedure is suggested. 

The time integral defined by Eq. (46) is evaluated 
(in general, by approximate methods) along the paths 
IABF,...,IKLF. Since each of these paths is asso- 
ciated with a different value of the constant //, the inte- 
gration process yields a result of the following type: 


ty = (81) 
The analogous integration for Eq. (45) leads to 
X, = (S82) 


Eq. (Sl) supplies the value of /7 (and therefore the 
special Eulerian arc) associated with the desired value 
for the time. In turn, Eq. (82) determines the maxi- 
mum (or minimum) distance which the rocket-powered 
aircraft can fly in the established time interval. 

(9.1.6) Numerical Example—With the object of illus- 
trating the previous theory, a numerical computation 
is carried out for the following set of conditions (see 
Fig. 7): 


V,/V, = 0.30 


V,/V. 


m, ms. = 0.074, 


0.38 


mz = 0.027, 


Due to the qualitative character of the example the 
simplifying assumption 6,,,; = © is used. The time 
t, and the distance X, are represented through the fol- 


lowing dimensionless parameters :* 
gt /EV., gX/EV? (83) 
E = 1/2VKCpo (84) 


where 


is the maximum lift-drag ratio of the aircraft. 


~ 
Ve 
| 
mV). 
0, X,=0 
at 


JOURNAL OF 


H=0 


H=-O1 


u,=030 , 4=038 


m, m 
m= 0074 , m= 0.027 


05 06 "9 O07 08 


Fic. 8. Maximum distance as a function of the overall flying 


time (parabolic polar with constant coefficients). 


w(m,V,H) > 


Fic. 9. Comparison between the Eulerian are /.1/NF and the 


arbitrary control path JAF. 


The solution of the boundary value problem is plotted 
(Fig. 7) in the parametric form embodied by Eqs. (S81) 
and (82). For a time interval ¢; = O.71EV,/g the 
trajectory maximizing the range is the path /R/F asso- 


ciated with the variable thrust sub-are /J = +0.1. 
Analogously, for a time interval ¢; = O.S7EI’, g the 
solution corresponds to the path /7 = —0.1.  Elimi- 


nation of the parameter // leads to an expression of the 
type X, = C(t,) which is indicated in Fig. 8. This 
graph directly supplies the maximum range achievable 
in a given time. 

From an engineering point of view it is of interest to 
remark that the maximum absolute range is obtained 
by using the sub-are /J = 0.** By flying along the 
sub-are JT = +0.12—-which corresponds to flight ve- 
locities about 20 per cent higher—-one may save about 
20 per cent in time with a penalty of only 5 per cent 
in range. Analogously, by flying along the sub-are 
IT = +0.20-—-which corresponds to flight velocities 
about 30 to 60 per cent higher than the speeds asso- 
ciated with 7 = 0--a saving of about 40 per cent in 
time results with a loss of 20 per cent in range. 

(9.1.7) Sufficiency Proof for the Solution—-Assume 
now that the special Eulerian path, satisfying the 
boundary conditions indicated in Paragraph (9.1.5), has 
been found. Assume also that 77 is the particular 
value of JZ which enables the aircraft to fly from J to 
F in the prescribed time ¢;. Two possibilities exist 
the case where /7 < 0 and the case where [7 > 0. 
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For [7 < 0 (Fig. 9) the proof of the maximum rang 
character of the trajectory /.1/.NF of Fig. 9 is carried 
out by comparing it with the control path JAF; on 
must exploit the properties of the function w an 
Green’s theorem as in Paragraph (9.1.2). In addition 
the inequality 0 < 8 < Byax must be kept in mind, be. 
cause it determines the class of admissible displace. 
ments for the aircraft in the (m, I’) plane.* 

An analogous technique can be used for the case 7 > 
0. One can prove that the Eulerian are /.\/.NF yields 
a maximum range solution if it is entirely enclosed jn 
the two regions A and B (Fig. 4); it yields, on the con- 
trary, a minimum range solution if it is entirely en- 
closed* in the two regions B and C. 


CONCLUSIONS 


The burning program which minimizes an arbitrary 
function of the final values of time, mass, distance, and 
velocity is investigated. It is shown that the totality 
of extremals is composed of zero thrust sub-arcs, sub- 
ares to be flown with maximum engine output, and 
variable thrust sub-arcs. 

Within the general frame of the previous approach, 
the thrust programing maximizing the range for a pre- 
scribed time of flight is discussed and a method for 
constructing extremal paths is supplied. 

By a suitable choice of the burning technique, the 
overall time can be considerably decreased with respect 
to the time associated with the solution of maximum 
absolute range. Numerical examples show that up to 
20 per cent savings in time can be achieved with a 
penalty of less than 5 per cent in range.f 
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* The writer has not yet found a satisfactory way for covering 
the more general case of ares which belong to the three regions 
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+ Extensive numerical computations have been carried out in 
reference 7 for the case where the coefficients of the parabolic 
polar depend on the Mach Number M. These computations 


have confirmed the above point of view 
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Experiments on Boundary-Laver Transition 
at Supersonic Speeds’ 


E. R. van DRIEST* anv J. CHRISTOPHER BOISON** 
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SUMMARY 


Tests were conducted in the 12-in. continuous supersonic wind 
tunnel of the Jet Propulsion Laboratory, California Institute of 
Technology, to determine the effects of surface cooling on 
boundary-layer transition at supersonic speeds. The effects of 
cooling were investigated at test section Mach Numbers of 1.97, 
281, and 3.84 with an internally cooled 20-in. 10° (apex angle) 
smooth cone in the presence of three levels of supply-stream tur- 
bulence (0.4, 2, and 9 per cent) and several single-element rough- 
nesses at fixed axial location. Transition data were obtained 
optically by means of a magnified-schlieren system. The results, 
for the range of Mach Number investigated, indicate that (1) 
transition on a smooth cone can definitely be delayed by surface 
cooling, (2) transition promoted by either supply-stream turbu- 
lence or surface roughness can also be delayed by surface cooling 
depending upon degree of turbulence or relative roughness re- 
spectively, and (3) the adverse effects of increased turbulence and 
roughness decrease with increasing Mach Number. 


SYMBOLS 


M5 = local Mach Number at outer edge of boundary 
layer 
T3 = local ambient temperature at outer edge of 


boundary layer 
wall temperature 


l../Ts = average wall-to-local-free-stream temperature ra- 
tio 

P = total pressure in supply chamber of wind tunnel! 

ps = fluid density at outer edge of boundary layer 

Ms = fluid viscosity at outer edge of boundary layer 

ug = velocity at outer edge of boundary layer 

u’ = root mean square velocity fluctuation 

(u'/ue) = intensity of turbulence in supply chamber of wind 
tunnel 

r = local recovery factor 

k = height of roughness 

Xx = distance from tip along cone generatrix 

Ny = distance from tip along cone generatrix to transi- 
tion 

x, = distance from tip along cone generatrix to rough- 
ness element 

6 = boundary-layer total thickness 

6, = boundary-layer total thickness at roughness ele- 
ment 

6* = boundary-layer displacement thickness 
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6,* = boundary-layer displacement thickness at rough- 
ness element 

5% = dimensionless displacement thickness 

Res = Reynolds Number in terms of conditions at outer 


edge of boundary layer and distance x 

Res, = Reynolds Number in terms of conditions at outer 
edge of boundary layer and distance to transi- 
tion x, 

Res, = Reynolds Number in terms of conditions at 
outer edge of boundary layer and distance to 
transition x7 for the smooth body at minimum 


turbulence 


Res, = Reynolds Number in terms of conditions at outer 
edge of boundary layer and height of roughness 
Re;, = Reynolds Number in terins of conditions at oute: 


edge of boundary layer and distance to rough- 
ness element x, 


Res/in. = Reynolds Number per inch at outer edge of 
boundary layer 

ims = insulated wall (zero heat transfer 

min. crit. = minimum critical 

turb = turbulence 


INTRODUCTION 


T WAS ORIGINALLY PREDICTED" * that, for supersonic 
flow over a smooth surface, a laminar boundary 
layer could be highly stabilized by drawing heat out of 
it at the surface of contact. Whether this phenomenon 
could be realized in free flight with practical surface 
roughness was open to question. For this reason, a 
fundamental experimental research program was initi- 
ated to investigate the control of transition by cooling 
in the presence of free-stream turbulence and surface 
roughness. The results would have practical import- 
ance for internal cooling systems in steady flight as 
well as naturally-induced cooling in transient flight. 


WALL-TO-LOCAL-FREE - STREAM 
TEMPEATURE RATIO, Tw/ Ts 


° ' 2 3 4 5 6 7 R 
LOCAL MACH NUMBER , Mg 


Fic. 1. Cooling required to stabilize a laminar boundary layer 
on a smooth cone. 
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Table | 
TABULATED RESULTS 
10° (Apex Angle) Cone 
M 3 (Re;/ inch) 
avg 
1.90 312R 0.54 X 10° 
Note: All trips are circular wires 
2.70 228 R 0.67 x 106 except as noted. 
3.65 152R 0.50 X 10° 
M Config Re 0-4 M Config te lee 0-4 Config. Ire “10-4 
8 8 Sr . 3 $ Ts br 
Turb. Trip Trip 
1.90 0.4% 1.628 5.841 1.90 0.002” 1.639 3.230 2.70 0.005” 1.853 4.833 
1.481 9.477 1.530 3.367 1.589 9.397 
1.459 10.270 1.394 3.258 2.249 4.087 
1.539 7.923 1.411 3.333 1.782 6.189 
1.494 9.110 1.636 3.267 1.660 8.186 
1,557 6.988 1.271 3.492 1.558 9.384 
1.518 8.438 1.003 3.357 1.513 11.979 
1.624 6.277 0.579 3.141 2.259 4.126 
1.521 8.471 1.420 3.540 1.807 6.698 
1.546 7.906 1.684 7.891 
1.595 6.702 Turb. 1.556 9.946 
1.518 9.176 2.70 0.4% 2.281 3.947 1.552 11.841 
1.486 8.987 1.816 6.069 Trip 
1.608 6.450 1.741 6.588 
1.617 6.062 1.797 5.618 — oan 2.263 4.356 
1.599 7.038 2.105 4.757 
1.876 5.437 1.972 4.809 
1.611 5.405 1.526 11.013 1.793 5.945 
1.586 6.062 1.775 6.736 1.748 6.546 
1.548 6.599 2.262 4.331 1.683 7.404 
1.521 7.065 2.116 4.500 1.589 7.50 
1.455 8.437 1.967 5.104 1.488 9.006 
1.483 8.010 1.814 5.428 1.504 9-38 
1.523 7.102 1.732 6.070 1.425 Pe 
1.569 6.286 1.687 7.335 pe 
1.490 7.938 1.593 7.988 oan canes 
1.455 8.124 1.528 9.464 
Lar 11.4688 1.259 10.063 
0.655 4.469 
2.261 4.316 
1.635 4.268 1.929 5.115 Trip 
1.370 5.124 1.762 6.528 2.70 0.002" 2.265 4.422 
1.195 5.266 1.680 7.316 2.112 4.555 
1.048 5.647 1.705 7.370 1.961 4.743 
0.669 6.119 1.626 8.149 1.797 5.611 
1.611 4.568 2.255 4.095 1.707 5.899 
1.517 4.676 1.708 7.282 1.682 6.750 
1.455 4.721 1,586 9.804 1.559 7.060 
1.377 5.189 1.538 11.873 1.595 7.888 
1.433 7.535 
Turb. 1.478 8.562 
1.622 6.203 2.70 2% 2.256 4.177 0.629 4.623 
1.608 6.396 2.081 4.508 1.191 7.370 
1.568 7.021 1.938 5.160 1.208 7.638 
1.551 7.565 1.809 5.822 1.950 4.668 
1.503 9.228 1.793 6.615 1.724 5.790 
1.491 9.496 1.729 6.930 1.492 7.014 
1.457 9.944 1.679 7.920 1.172 6.469 
1.637 8.070 0.969 6.370 
1.605 8.382 0.890 6.045 
1.641 4.724 1.540 11.220 
1,546 4.724 Trip 
1.468 4.724 Turb. 2.70 0.004" 2.321 3.996 
1.389 4.724 2.70 9% 2.254 3.349 2.190 3.960 
1.390 4.724 1.953 4.162 2.064 3.893 
1.318 4.724 1.615 5.469 1.904 3.893 
1.648 4.616 1.462 5.760 1.758 4.161 
1.325 4.675 1.786 4.776 1.614 4.087 
1.288 4.717 1.483 4.221 
1.133 4.474 Trip 0.825 3.551 
1.185 4.648 2.70 0.005” 2.264 4.384 0.808 3.624 
1.633 4.676 2.136 4.317 0.795 3.625 
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Table 


TABULATED RESULTS 
20’”— 10° (Apex Angle) Cone 


Ms Ts (Re; /inch) 
avg ovg 
1.90 312R 0.54 X 10° 
Note: All trips are circular wires 
2.70 228 R 0.67 X 10° except os noted. 
3.65 152R 0.50 X 106 
Confi Tw Ike “10-47 Config. ke. “10-4 M Config Re. *10-4 
Ms g- T | 8 T5 3 8 
Turb. Trip 
3.65 0.4% 3.352 3.245 3.65 0.004’ 2.107 5.788 
2.687 4.515 1.925 7.435 
2.289 5.004 
1.823 7.484 Trip 
2.008 6.268 3.65 0.008" 3.321 3.337 
3.307 3.387 2.952 3.838 
3.363 3.718 2.544 4.303 
2.967 3.855 2.351 4.739 
2.445 4.911 2.267 4.849 
2.077 5.977 2.113 5.385 
1.790 8.302 2.004 5.988 
3.292 3.277 
2.439 4.757 Trip 
2.043 6.272 3.65 0.010" 3.313 3.419 
2.252 5.061 2.972 3.629 
1.907 6.748 2.691 5.934 
2.770 3.986 2.467 4.444 
2.213 4.840 
Turb. 2.097 5.093 
3.65 9% 3.294 3.312 1.938 5.370 
3.316 3.498 1.849 6.053 
3.065 3.571 1.672 7.978 
2.415 4.541 2.847 3.749 
2.016 6.020 2.624 4.224 
1.739 7.512 2.240 4.709 
Sharp 
Trip Trip 
3.65 0.001" 3.298 3.379 3.65 0.0105" 3.301 3.294 
3.176 3.488 1.961 5.412 
2.898 3.923 
2.643 4.426 Trip 
2.774 3.812 3.65 0.012" 3.361 2.396 
2.947 3.929 3 Di- 3.196 2.347 
3.094 3.754 mension 3.041 2.513 
2.488 4.891 2.819 2.594 
2.351 5.289 2.893 2.336 
2.234 5.422 1.936 2.515 
2.060 6.020 1.226 2.300 
1.985 7.386 
Trip 
Trip 3.65 0.020’ 3.336 3.041 
3.65 0.002"’ 3.323 3.532 2.864 3.250 
3.202 3.731 2.457 3.407 
3.105 3.731 2.111 3.822 
2.929 3.930 1.540 4.133 
2.772 4.267 0.998 4.107 
2.597 4.676 0.898 3.556 
2.499 4.962 
2.362 5.061 Trip 
2.189 5.316 3.65 0.037" 3.357 2.552 
2.027 6.557 3.054 2.558 
1.975 7.508 2.580 2.660 
1.926 2.721 
Trip 1.030 2.396 
3.65 0.004’ 3.417 3.629 
3.220 3.640 
2.922 3.904 
2.658 4.341 
2.525 4.875 
2.404 5.024 
2.337 5.672 
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4.833 a 
9.397 
4.087 
6.189 
8.186 
9.384 
11.979 : 
4.126 
6.698 
7.89) a 
9.946 
11.841 
4.356 
4.757 a 
4.809 
5.945 
6.546 
7.604 
7.517 
9.006 ‘ 
9.338 
9.870 
10.005 
10.405 
10.005 
10.063 
4.469 
4.422 
4.555 
4.743 
5.611 
5.899 
6.750 
7.060 4 
7.888 
7.535 
8.562 | 
4.623 
7.370 
7.638 
4.668 
5.790 
7.014 
5.469 
5.370 
5.045 
}.996 
960 
4 
893 
161 
.087 
221 
624 
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2. View of 20-in. 10° cone in supersonic wind tunnel. 
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Fic. 3. Schematic of instrumentation and cooling system. 


The cooling required to stabilize the laminar bound- 
ary layer on a cone in an unheated wind tunnel to 
various minimum critical Reynolds Numbers is shown in 
Fig. 1. 
cal Reynolds Numbers are but estimates, the infinite- 
stability curve is quite accurate. Particularly signifi- 
cant in the Figure is the observation that over the range 
of Mach Number of about | to S it is possible to stabilize 
the boundary layer completely with a practical cooling 
Also significant is the observation that the stabil- 


Although the lower values of the minimum criti- 


rate. 
ity of the boundary layer is very sensitive to tempera- 
ture change as the complete-stability curve is ap- 
proached. The complete-stability curve chosen as a 
base from reference 2 represents temperature conditions 
in an unheated supersonic wind tunnel where it might 
be expected that the temperature within the boundary 
layer would be such that uw « Tand Pr ~ 0.75. In free 
flight, higher base curves would be used, owing to the 
change from « tow T' as the temperatures 
throughout the boundary layer are increased. 

This report describes experiments conducted at test- 
section Mach Numbers of 1.97, 2.81, and 3.84. The 
body tested was a 20-in. 10° (apex angle) cone. In 
order to obtain wall-to-local-free-stream temperature 
ratios required for strong stabilization, the cone was 
cooled internally with premixed liquid and gaseous 
nitrogen. The range of cooling is indicated in Fig. 1. 
Appropriate screens were installed in the supply cham- 
ber to produce turbulence intensities there of 0.4, 2, and 
9 per cent. Various sizes of circular wires were used as 
single-element roughnesses located 3 in. from the nose of 


1957 


DECEMBER, 


the cone. While the circular wires comprised the by 
of the roughnesses, one three-dimensional trip was j 
vestigated. Transition was located optically by m 
nifying the schlieren image normal to the bound; 
layer. 

The principal result of these experiments was thy 


ay 


cooling was definitely effective in delaying transitio; 
depending upon the intensity of turbulence and th 
magnitude of the roughness. Another important resy) 
was that the transition- promoting effects of turbulen 
and roughness diminished with increase in Mach Nup 
ber. The single three-dimensional trip introduced wa 
as effective in promoting transition with cooling as, 
circular trip three times as high. The data wey 
studied in the light of the recent correlation work 6 
Dryden for incompressible flow.’ 

Table | gives all data on transition presented in thi: 
paper. 


APPARATUS 


Wind Tunnel 


The experiments were performed in the 12-in. super- 
sonic wind tunnel of the Jet Propulsion Laboratory 
(CIT). The tunnel! is of the open circuit, variable pres- 
sure, continuous operation type. A 
Mach Number from approximately 1 to 4 is made 


variation of 


available through the use of a flexible-wall nozzle. Be- 
tween Mach Number 1.0 and 2.2 the test section is 12 
in. square; above Mach 2.2 the test section height is re- 


duced to 9in. During the tests the supply temperature 
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Fic. 4. Typical schlieren photos showing the effect of several 


cooling rates for the smooth 10° cone. \/5 = 3.65. 
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Typical surface-temperature distributions for several 


cooling rates on a smooth 10° cone. 6 = 3.65. 
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distribution along a smooth 10° cone without heat 
Ms = 1.90. 


Schlieren photo of transition compared with tem- 
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Fic. 7. Schlieren photo of transition compared with tem- 


perature distribution along a smooth 10° cone without heat 
transfer. 1/5 = 2.70 
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Fic. 8. Schlieren photo of transition compared with tem- 
perature distribution along a smooth 10° cone without heat 
transfer. M5 = 3.65. 
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varied from 70°F. to 100°F. with a dew point variation with 1 S-in. wall thickness. 321 stainless was chose} ~_ 
of —10°F. to —40°F. A Reynolds Number per inch of for three reasons —first, it does not precipitate j, 
approximately 0.5 X 10° was obtained on the 10° carbon when exposed to high welding temperature 
conical test body for the three Mach Numbers investi- second, it is noncorrosive, and third, it has desirabj, 
gated. thermoelectric properties. The model was fabricate 
The turbulence level in the tunnel was varied by in- in two halves split along the horizontal plane main] 
troducing various grid and screen combinations into to allow installation of inner-surface thermocouple 
the supply-chamber section. The normal-operating tur- All thermocouples were installed at the same time. Th 
bulence level of the tunnel was about 0.4 per cent (re- two sections were then welded together and the surfag 
ferred to in this report as minimum turbulence). A ground and hand polished to a nominal surface finish ¢) 
ribbon-mesh combination was used to raise the turbu- approximately 10 microinches. 
lence level to about 9 per cent (maximum turbulence). A spewing or distribution head was made from severg 
These two levels were determined experimentally by sections of l-in. to 1 S-in. brass tubing arranged in , 
Dr. J. Laufer of JPL. It was possible, therefore, to stepwise fashion. This head was drilled with a larg 
change the supply-turbulence level by a factor of ap- number of orifices ranging from 0.018 to 0.092 in. in 
proximately 20. <A third level of turbulence (2 per cent) diameter, attached to the sting support and inserted 
investigated was determined by empirical methods. into the outer conical shell assembly, thus providing the 
Test Bod mechanism for coolant distribution. 
we . The conical shell was instrumented with a total of 52 
The test body of the present investigation was a 20- thermocouples, the location of which were as follows 
in. 10° (apex angle) conical shell of 321 stainless steel 39 on the outside-surface upper generatrix, spaced 0.5 
12 
10 Frc 
fo) 6 
= Res/INCH = 0.54-10 
| | 
| 
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Ms =3.65 | | 
x | Res /INCH~ 0.50-10° | | | 
| | | 
= | | 
5 | Ms =2.70 | | | 
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| | | } 
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S 4 
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| 3.65 | | 
le | 
1.0 1.4 1.8 2.2 2.6 3.0 3.4 38 
AVERAGE WALL~TO-LOCAL~STREAM TEMPERATURE RATIO, Tw/Ts 
Fic. 9. Effect of surface cooling on transition Reynolds Number for several local Mach Numbers on a smooth 10° cone F 
(0.4 per cent supply turbulence). 
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in. apart with the first thermocouple located 3.4 jp, 
from the tip; 12 on the inside-surface upper generatrix 
spaced 1.5 in. apart with the first thermocouple located 
3.7 in. from the tip; and 5 on the outside-surface lower 
generatrix, spaced 4-in. apart with the first thermo. 
couple located 3.8 in. from the tip. The latter fiye 
thermocouples were used for checking the circumferen. 
tial temperature symmetry. 

For the surface-roughness investigation, a series oj 
circular wires ranging in diameter from 0.0005 to 0.037 
in. was used as single-element boundary-layer trips 
The trips were accurately constructed to fit the conical 
shell at 3 in. from the tip. 

Fig. 2 shows the model in place in the wind tunnel, 
A wire trip can be seen near the tip. 


Cooling System 


The cooling system producing the low temperatures 
necessary for these tests utilized both liquid (LN.) and 
gaseous (GN,) nitrogen. By mixing these 
ponents it was possible to maintain steady-state tem- 


com- 


peratures from room temperature down to —320°F. on 
the model. 

Fig. 3 gives a schematic diagram of the cooling sys- 
tem. The liquid nitrogen was stored in a stainless steel 
sphere capable of pressurization up to 400 psig. The 
gaseous nitrogen was supplied from a manifold oi 4s 
bottles at 2,250 psig. Prior to entry into the model, the 
liquid was passed through a spray nozzle and mixed 
with GN,. The cooled fluid entering through the sting 
mount was distributed inside the conical shell by a dis- 
tribution head and finally exhausted through the model 
base into the wind-tunnel stream. 

The sphere was mounted on a platform scale so that the 
amount of LN, in the container could be determined and 
the approximate flow rates calculated. LN, flow rates 
as high as 2 gal. per min. were used in the present tests. 

Model surface-temperature control was obtained 
through regulation of both the LN, and GN, supply 
pressures. LN, and GN, pressures were as high as 
approximately 200 psig. 


Thermocouple Instrumentation 


A schematic of the instrumentation used for tempera- 
ture measurement and recording is shown in Fig. 3. 
Model temperatures were measured at 52 locations, 
using a common lead system. The cone and a single 
321 stainless steel lead served as the common leg while 
constantan leads, silver-soldered to the various cone lo- 
cations, composed the other legs. All the constantan 
leads passed through two channels between the spewer 
head and the conical shell (immersed in the coolant to 
minimize conduction in the constantan leads). 


The recording of temperatures was divided between * 


two Leeds and Northrup Speedomax recorders using 
identical circuitry for each instrument. The leads 
passed from the model through a thermocouple push- 
button-type selector switch, a pulsing mechanism, and 
a Brown potentiometer, into the recording instrument. 
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sing mechanism was used to facilitate identifica- 
the successive thermocouple outputs. The 
meter was used to introduce a bucking voltage 


The pul 
tion of 


potenti 
to position the recorder zero-points advantageously, 
thus making it possible to record the complete tempera- 


ture range investigated in the experiments. A Dewar 
flask filled with cracked ice and water was used for 


maintaining a 39°F, reference junction. 
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It was possible to record all the thermocouple emf’s 
for each run within a period of about 2 min. The ac- 
curacy of temperature measurement was within +1°F. 

The constantan leads for measuring outer surface 
temperature were directed from the inside through small 
holes drilled in the shell. A portion of the hole was 
packed with a ceramic type of insulation and the lead 
then silver-soldered from the outside; this insured a 
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good surface temperature reading. The inner surface 
thermocouples were silver-soldered in a small conical 
indentation of the shell’s inner wall. 

The thermocouples were calibrated at the beginning 
of the test schedule by placing the model in a very 
accurately temperature-controlled “cold box’ and the 
outputs compared to the reading obtained on several 
standard thermocouples located near the model. 
Thermocouple outputs were measured on a K-2 Leeds 
and Northrup potentiometer, 


Magnified Schlieren System 

A special optical system? utilizing a cylindrical lens 
was constructed for the purpose of locating transition. 
This system was placed just aft of an auxiliary knife 
edge of the JPL schlieren system and aligned parallel to 
the upper generatrix of the cone surface. The lens 
magnified the distance normal to the cone surface by a 
factor of about 19, thus distorting the shape of the 
boundary layer and allowing the location of transition 
to be determined accurately. A 16-mm. movie camera 
was used to photograph the boundary layer as seen 
through the special lens system. The JPL standard 
schlieren system was utilized only for the purpose of 
obtaining front-lighted schlieren photographs. The 
front-lighting technique was especially desirable in order 
to record and evaluate the extent of frost formed on the 
model during excessive cooling (to be discussed later). 
During the period of the recording of temperature data, 
it was possible to take both magnified as well as front- 
lighted schlieren photographs. 


TEST PROCEDURES 


The test schedule was carried out in two main 
phases, namely—phase 1, the effect of cooling on 
boundary-layer transition for a smooth cone with 
several levels of free-stream turbulence (0.4, 2, and 9 
per cent), and phase 2, the effect of cooling on bound- 
ary-layer transition under the influence of specific 
ordered surface roughnesses (single element) with 
minimum free-stream turbulence. 

The above tests were conducted at the highest possi- 
ble Reynolds Number per inch in order to achieve 
transition in a most forward position on the model. 
This procedure allowed the observation of relatively 
large movements of transition by cooling. The two 
phases were investigated at test-section Mach Numbers 
of 1.97, 2.81, and 3.84 (corresponding to local Mach 
Numbers .1/; at the outer edge of the boundary layer of 
1.90, 2.70, and 3.65, respectively) with nominal local 
Reynolds Numbers per inch of 0.54 X 10°, 0.67 X 10°, 
and 0.50 X 10°, respectively. 

The general testing procedure first reached a steady- 
State zero-heat-transfer condition, whereupon tempera- 
ture data and magnified, as well as front-lighted, 
schlieren photographs were taken. The model was then 
cooled in decrements. Owing to the sensitivity of 
transition to cooling, small temperature decrements 
were necessary. For example, decrements as small as 
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5°F. were used in the case of 1/; = 1.90. At each pre- 
determined temperature, steady-state conditions were 
obtained before data were taken. The temperature of 
the model was lowered until transition was either 
moved off the model or maximum cooling had been 
reached. 

A separate series of experiments with zero heat 
transfer were carried out to observe the effect of 
Reynolds Number per unit length on transition. 

The range of cooling achieved during the conduct of 
the investigations is shown in Fig. 1. It is seen that 
local temperature ratios could be attained within the 
area of complete stability for all Mach Numbers in- 


vestigated. 
RESULTS AND DISCUSSION 


Location of Transition 


As mentioned above, the location of transition was 
determined optically in these experiments by means of 
a specially-designed magnified (normal to the flow) 
schlieren apparatus. Typical photographs used for 
measurement are shown in Fig. 4+ for local Mach Num- 
ber 3.65 and various rates of cooling. The surface- 
temperature distributions corresponding to the photos 
are given in Fig. 5. It is seen that transition can be im- 
mediately identified in the photos, while the location 
from surface-temperature distribution alone is inde- 
terminate except for the zero heat-transfer condition. 

Originally, it had been intended to determine transi- 
tion from heat-transfer-coefficient calculations, and for 
that reason thermocouples were installed at the inner 
surface as well as at the outer surface of the model. 
However, this method was soon abandoned in favor of 
the magnified schlieren system. 

It is interesting to compare some schlieren photos for 
zero heat transfer and the corresponding temperature 
distributions to find out what part of the temperature 
distributions corresponds to the location of transition 
as seen in the photos. Figs. 6, 7, and 8 show typical 
schlieren photographs and corresponding temperature 
distributions for the zero heat-transfer case. From these 
Figures it is apparent that the photographic transition 
point correlates better with the top than with the bot- 
tom of the temperature rise, and this correlation was 
found to be true for all zero-heat-transfer data. 


Transition on a Smooth Body 


Transition data for a smooth cone in a wind tunnel 
with minimum turbulence are plotted in Fig. 9 as a 
function of Mach Number and cooling rate. It is evi- 
dent from the Figure that cooling delays transition and 
that, for zero heat transfer, transition decreases with 
increase in Mach Number for the range of Mach Num- 
bers investigated. With increased cooling the curves 


seem to approach approximately the lines indicating 
complete stability for infinitestimal two-dimensional 
disturbances. The three curves of Fig. 9 will be the 
bases relative to which the effects of turbulence and 
roughness will be measured. 
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The average wall-to-local-free-stream temperatur 
ratio 7,,/7; was the weighted average over the area uy 
to the point of transition as observed from the photo. 
graphs. 


Effect of Supply-Stream Turbulence on Transition 


Figs. 10, 11, and 12 show the effect of supply-stream 
turbulence on transition with cooling for a smooth cone 
These Figures indicate clearly how an increase in Mach 
Number diminishes the transition-promoting effect oj 
turbulence; thus the strong effect of a 9 per cent sup. 
ply-turbulence at local Mach Number 1.90 has prae 
tically disappeared at Mach Number 3.65. 

How transition varies with Mach Number and in- 
tensity of supply turbulence for zero heat transfer js 
seen in Fig. 13. 

Crossplots of Figs. 10-12 are shown in Figs. 14-17 in 
order to indicate separately the effects of Mach Number 
and cooling relative to the minimum-turbulence smooth- 


body data. 


Effect of a Single Roughness Element on Transition 


Single wire rings of various sizes were placed at a 
position of 3 in. from the nose of the cone. While most 
wires were circular cylinders, a sharp-edged and a three- 
dimensional type were tested. 

Fig. 1S shows an ordinary spark shadowgraph of the 
flow past a circular wire at a Mach Number of 3.65. 
Undoubtedly the flow pattern associated with the 
roughness at supersonic speeds strongly influences the 
position of transition. 

Figs. 19-21 give the transition data for the cooled 
cone with single-element surface roughness. As with 
supply-stream turbulence, the effect of increase in 
Mach Number is to diminish the effect of the roughness 
disturbances. Also, cooling is effective in delaying 
transition, depending upon the size of the wire. The 
strong reversal in transition with increased cooling as 
seen in Fig. 20 for the 0.001-in. and 0.002-in. trips can be 
accounted for through the argument that, as cooling 
proceeds, the stabilizing effect of cooling on the bound- 
ary layer, fore and aft of the trip, is eventually offset 
by the increasing disturbance of the trip, owing to the 
increase of the size of the trip relative to the boundary- 
layer thickness. This phenomenon should also occur at 
M,; = 1.90 and 3.65 for intermediate trips. 

Figs. 22-25 are crossplots of Figs. 19-21 with the 
ratio of the roughness height to boundary-layer dis- 
placement thickness at the roughness position, k 6,*, as 
an independent variable, following the method of corre- 
lation suggested by Dryden for incompressible flow.’ 
These curves show the separate effects of Mach Number 
(Fig. 22) and cooling (Figs. 23-25) on transition relative 
to the smooth body data with minimum turbulence. 
The spectacular role of Mach Number in damping the 
effect of roughness on transition is evident in Fig. 22, 
where it is seen that, with increasing Mach Number, 
increasingly large ratios of roughness height to bound- 
ary-layer displacement thickness were necessary to 
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Fic. 30. Frost patterns as recorded by front-lighted schlieren. 
promote transition relative to the smooth-body condi- 
tion. (Note that for 1/; = 3.65 and insulated wall 
a roughness-height-to-boundary-layer-displacement- 
thickness ratio of 4 corresponds to a roughness-height- 
to-boundary-layer-total-thickness ratio of 2! ,). 

It is also important to observe from Fig. 21 that the 
three-dimensional roughness of size 0.012 in. promoted 
transition to the same degree as a circular wire of 0.037 
in. 

The dimensionless displacement thickness for a lami- 
nar boundary layer on a cone as a function of Mach 
Number and heat transfer is given in Fig. 32 as com- 
puted by means of the Crocco method in reference 6. 


Effect of ‘‘Reynolds Number per Unit Length’”’ 


For a smooth body, if the Mach Number, temperature 
ratio, and intensity of turbulence are fixed, it is ex- 
pected that the transition Reynolds Number 
would remain fixed, and, therefore, x7 would vary in- 
versely with psuv5/u;. In order to prove this expecta- 
tion, a separate series of experiments were conducted in 
which at each Mach Number the Reynolds Number was 
varied by changing the supply pressure—i.e., density. 
That the actual distance to transition follows the inverse 
law quite well is evident in Figs. 26, 27, and 28. How- 
ever, such deviations as in Figs. 27 and 2S have not been 
fully explained. Fig. 29 shows transition photos corre- 
sponding to some of the data points in Fig. 26. 


1957 


Frost Formation 


A frost formation appeared on the model when th 
model was sufficiently cooled. Since there may be some 
question about the effect of the frost formation on the 
data obtained, a detailed dicussion of the general phe. 
nomenon is in order. 

As the model was cooled, frost first began to appear at 
decreasingly lower temperatures of approximately 
—SO°F. to —115°F. with increase of local Mach Num 
ber from 1.90 to 3.65. The first traces of frost became 
evident in the form of a fine ring at approximately 6 in 
from the tip. With increased cooling, the frost pattern 
would spread until it covered practically the entire 
model, except for the first 3° , inches. Typical tem- 
perature distributions for local Mach Number 3.65 are 
shown in Fig. 5. In that Figure it is seen that the 
model reached its initial temperature of —115°F. (cor. 
responding to 7,,/7'5 = 2.30) at about 6 in. from the tip 
of the cone and that at the greatest cooling rate the 
frosting temperature of —115°F. moved to 4.2 in. from 
the tip. It is further seen that, owing to the strong 
temperature gradient near the nose, it was quite im- 
possible to produce frost within about 3° , in. from the 
tip. This nose temperature gradient occurred as a re- 
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perature distributions of Figs. 6, 7, and 8. 
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sult of the larger heat-transfer coefficients of the 
boundary layer and the inadequacy of the cooling sys- 
tem within the model near the tip to carry the heat 
away. It is important to emphasize that frost never 
covered the trips, since the trips were located 3 in. from 
the tip and remained, therefore, in a region devoid of 
frost. 

The possible effect of the frost formation on the ex- 
perimental results will now be discussed. Fig. 30 
shows front-lighted schlieren photographs which indi- 
cate three circumstances under which transition could 
be affected. The three conditions, A, B, and C, corre- 
spond to data points marked in Figs. 20 and 21.  Pic- 
ture A shows sufficiently frost on a smooth cone that a 
shock wave is formed at the leading edge of the frost. 
However, according to the position of the corresponding 
data-point A in Fig. 20, it appears that the frost front 
was not of sufficient height to cause an appreciable 
change in the trend of the smooth-body curve of Fig. 20. 
Picture B shows heavy frost on a smooth cone, but with 
no leading-edge shock. Since the corresponding data- 
point B in Fig. 21 also seems to fall on what may be 
called a smooth-body curve, apparently the crystalline 
structure of the frost was so fine that it did not con- 
tribute to the promotion of transition. Picture C 
shows frost aft of the 0.010-in. trip without a frost 
shock. Again, the corresponding data-point in Fig. 21 
indicates no abnormal deviation from the established 
trend of the 0.010-trip data as cooling proceeds. 

From the above observations and remarks, it seems 
safe to conclude that the frost did not interfere appre- 
ciably with transition. At any rate the presence of the 
frost problem is another argument for heated wind 
tunnels or the use of working fluids composed of com- 
ponents which remain gaseous at low temperatures. 


Recovery Factor 


Since the recovery factor for insulated surfaces can 
be readily computed, given the temperature distribu- 
tion and the Mach Number, and since some workers are 
interested in these values, such calculations were car- 
ned out for the temperature distributions of Figs. 6, 7, 
and 8 and presented in Fig. 31 for completeness. 


CONCLUSIONS 
Supply-stream turbulence and single element surface 
roughnesses are decreasingly effective in promoting 
transition as the Mach Number is increased. 
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Displacement thickness on a cone as a function of free 
stream Mach Number and heat transfer 


Fic. 32. 


Cooling is effective in delaying transition regardless 
of whether transition is promoted by turbulence or 
roughness, depending upon Mach Number, intensity of 
turbulence, and height of roughness relative to the 
boundary-layer thickness. 
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SUMMARY 


An approximate method of solution for the supersonic flow 
about axisymmetric bodies is given. It is based on the method 
of linearized characteristics and is shown to provide a rapid and 
The method has the advantage 
As a nu- 


accurate means of calculation. 
that its limits of applicability are easily determined. 
merical example, the flow about the nose of the RM-10 body has 
been analyzed for several flight Mach Numbers and the results 
shown to be in good agreement with other methods of analysis. 


SYMBOLS 


= sound velocity 

= unit perturbations defined by Eqs. (16) 

= initial ratio of shock to body curvature 

= parameters in the shock equation 

= nondimensional Cartesian coordinates referred to 

the body length 

points in the flow field 

= Mach Number 

= nondimensional radial distance referred to the body 
length 

= entropy 

= intensity of the velocity vector 

= shock inclination 

= ratio of specific heats 

= parameters in the body equation 

= Mach angle 

= local inclination of the characteristic lines to radial 
lines (see Fig. 4) 

= inclination of the velocity vector to the x-axis 

= polar coordinate 


Subscripts 


b = quantities computed at the surface of the actual 
body 

= quantities computed at the surface of the basic 
cone 

properties of the superposed linearized perturbation 
fields 

= quantities referring to the basic flow field 


INTRODUCTION 


4 I ViiE PRACTICAL USEFULNESS of approximation meth- 
ods of solution for given problems is obviously 


Received April 18, 1957. 

7 The research was supported by the United States Air Force 
through the Office of Scientific Research of the Air Research and 
Development Command, under Contract AF 18(600)-693, 
Project No. R-352-30-12. Part of the material included in this 
report, including the basic ideas on the determination of the 
limits of applicability, were presented first in March, 1955.! 

* Research Group Leader. 

** Professor of Aerodynamics and Director of Aerodynamics 
Laboratory. 


900 


The Axisvmmetric Supersonic Flow Near the 
Nose of a Pointed Body of Revolution’ 


LUIGI G. NAPOLITANO* ann ANTONIO FERRI** 
Polytechnic Institute of Brooklyn 


indicated by their rapidity and accuracy and by the 
possibility of determining their limits of applicability 
Several works have been devoted to the solution of axi- 
symmetric flow fields around bodies of revolution.?~ 
However, lengthy numerical calculations are required 
if higher-order terms in the axial coordinate x are con- 
sidered. Furthermore, a certain amount ditii- 
culty in dealing with an encountered singularity seems 
to be present,'! whereas not enough concern has ever 
been given to a rapid and feasible way of determining 
the limits of applicability of the respective methods. 
The approach suggested here is believed to be rapid and 
accurate and readily lends itself to the determination of 
its limits of applicability. 


of 


The method is based on the linearized characteristics 
approach given in reference 7. The supersonic flow 
about a body of revolution at zero angle of attack is 
represented by a basic nonlinear flow field plus small 
perturbation fields which are linearized. The pertur- 
bations will propagate along the characteristic net oi 
the basic flow field’ and can be readily determined by 
finite difference calculations once their values at ap- 
propriate points are determined from the boundary 
conditions which cause the complete flow to differ from 


the basic flow. 


In the present application the basic flow has been 
taken to be that corresponding to the flow about a cone 
coincident with or close to that tangent to the body 
at its apex; thus, for many practical bodies, the tabu- 
lated values of reference 8 for the supersonic flow about 
cones can be used. The body shape is represented in 
terms of a power series in the axial coordinate, and the 
entire flow field is determined in terms of a similar 
series. The boundary conditions are applied at the 
surface of the basic shock and of the basic body. As 
was shown in reference 7 for problems wherein the basic 
flow is conical, the usual finite difference calculations 
required for the determination of the perturbations 
represented by any term of the series can be carried out 
in a thin strip bounded by a characteristic line of one 
family and by segments of the other family character- 
istic lines. The advantage of this fact with respect to 
the rapidity of the method is self-evident. 


It is of interest to note that once the unitary pertur- 
bations are determined for a given basic flow, they can 
be applied to the determination of the flow over a series 
of different bodies, provided that their associated basic 
flow field is the same. 
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OUTLINE OF THE METHOD 


Let the shape of the body in meridian planes be 


known and expressible as 


= Nov oa x? +1 (1) 
j=l 
with the hypothesis that \,> < \, for every Aj. The 
shock generated by this body is assumed to be similarly 
expressible in the form (see Fig. 1) 


n 
=kx + kw?! (2) 

j=0 
where kx is the basic shock generated by the basic cone, 
tangent or close to the one tangent to the body at the 
apex. Eq. (2) implies an attached shock which in turn 
requires that the flow be everywhere supersonic. The 
constant parameters k; are unknown and, as a measure 
of the deviation of the actual shock from the basic one, 
are related to the derivatives of the physical properties 
of the flow field at the nose. The flow field within the 
boundaries given by Eqs. (1) and (2) must be consistent 
with a uniform free stream in front of the shock given 
by Eq. (2) and will be accordingly described by 


ll 
| 


g=¢ < (3) 


Barred quantities are relative to the basic conical 
flow field as defined previously, and I’,, ¢), S;, are com- 
ponents of superposed perturbation flows. Consistent 
with the assumptions concerning the magnitude of the 
shock parameters k,, the squares of these components 
will be considered negligibly small. The quantities S; 
are uniquely related to the velocity fields since, to a 
linear approximation 
— 1)] grad S, = V-V, grad S$ + 

The gradients of S,; along and in the direction normal to 


the streamlines of the actual flow can be computed 
with respect to the streamlines of the basic flow.’ If 


901 


the basic flow is irrotational, the quantities S, will be 
constant along the streamlines of the basic flow. 

It has been shown’ that small disturbances propagate 
along characteristic lines of the basic flow; therefore, 
if I’, and yg; have been determined at any two points 
II and III of the basic characteristic net, the corre- 
sponding flow quantities at a point I, also a point of the 
basic net, can be found from the following equations :’ 


(1 + + — ZL) (¢)1 + a(S) = 
(1 — (Vidar (m L) (¢ > 


(+r) + (b + = 
(1 + (p + q) (¢ yu + 


where entropy terms are to be determined as functions 
of \,’'s from the projections of Eq. (4) along and in 
direction normal to the streamlines of the basic flow. 

The coefficients of Eqs. (4) and (5) are functions only 
of the known basic flow variables T and ¢ and their 
explicit expressions are given in reference 7. The 
initial values of 1’, and ¢, necessary to start the calcula- 
tions are found by relating the perturbation velocities 
and flow deflections along the basic shock to the un- 
known parameters .. 

Consider Fig. 1. The total velocity I’, at a point on 
the actual shock can be expressed as 


= Vp (dV dy) p, (yp (6) 


where P,; is on the basic shock and the conical char- 
acter of the basic flow field has been taken into account. 
Within the accepted linear approximation and by 
introducing Eq. (2) to define the quantity (¥, — ¥), 
Eq. (6) becomes 

n 
Vp = (T)p, (Vip, + 


(dV (1 + (7) 


The quantities (T7)»p, and (dT dy)p, are constant along 
the line Y = J and are given by conditions behind the 
basic shock with angle Y. These quantities will be 
indicated by the subscript os. The velocity I’p must 
also satisfy the shock relations corresponding to the 
perturbed shock ; thus, 


Vp = + (07 /06),, X 
tan (j+ | — tan (S) 


where (61 de),, is the derivative of the velocity behind 
the shock with respect to the angle of the shock evalu- 
ated ate = tan-'k = ¥. To the linear approxima- 
tion, Eq. (S) becomes 


Vp = V,, + (OF X 
| G+ (1 + | (9) 


j=0 
and equating Eqs. (7) and (9) yields the desired rela- 
tions as 
= (1 + + 1) 
(OV /de) — (dV (10) 


{| a 
#4 
-~ACTUAL BODY y=3);x’ 
: 
4 
yy 
Fic. 1. Geometry of the problem. (5) : + 
| n 
| a 
j=0 
| j=0 
S=S+)3> S?<S, 
j=0 
| 


902 JOURNAL OF THE AERONAUTICAL SCIENCES—DE 


since (V’)p, = (Vis)p, and xp = xp, + O(R,). 
Initial values for ¢; are similarly obtained as 


(¢)/Rw)p, = (1 + [G+ 1) X 
(O0¢/0€) — (11) 


As may be expected, the quantities j/k,’ and 
are constant along the basic shock. 

These initial values are all independent of the radial 
coordinate R. This reduces considerably the calcula- 
tions required to solve the flow field because both j/k jx’ 
and ¢, kx’ are only functions of y throughout the field. 
Then in order to obtain a solution of the problem, it is 
sufficient to choose the points of the basic character- 
istics net along radial lines with an increment in y 
selected according to the desired accuracy. Once the 
perturbation quantities have been determined at point 
B; (see Fig. 2) from any two points B; and By on the 
basic shock, the same values apply at point By, on the 
intersection of the radial line through B,; and of the 
basic characteristic line through B;. Thus, the re- 
quired elements are available to compute the pertur- 
bation flow properties at B; and, consequently, through- 
out the field until the radial line corresponding to the 
surface of the basic cone y = y, is reached. 

To complete the solution of the problem, the co- 
efficients k; must be determined. This is accomplished 
by satisfying the boundary condition that the actual 
body, within the linear approximation, be a stream- 
line. The flow deflection along the basic body is given 
by 


n 
go = (Se + (¢, (12) 
ju 
where the subscript ¢ signifies values computed at the 
surface of the basic cone. It might be pointed out in 
connection with Eq. (12) that the derivatives of the 
basic flow are finite at the surface of the cone. 
To the accepted order of approximation, one can 
write 


tan g = (tan ¢), + 


dey (1 + tan? g), (13) 


j=0 


and this expression must be equal to the corresponding 


one derived from the body Eq. (1). Accordingly, 
FAMILY 
—---SECOND FAMILY 
y=CONSTANT 


Fic. 2. Network for characteristic calculations. 
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Initial ratio of shock to body curvature against the 
combined supersonic-hypersonic parameter. 


Fic. 3. 


wt G+ = 


(tan 2). + (¢)/kix’) Ry’ [1 + tan? 


j=0 
which is satisfied for any value of x if 
ky = [Xo — (tan [1 + tan? =], 
ky = + [1 + tan® 2]. 


(14) 


When the basic cone is tangent to the body at its 
apex, (tan ¢), = Ay and &p is identically zero. 

When the parameters k; are known, the entire flow 
field is known. It is pointed out that the general prob- 
lem is thus completely solved by means of simple finite 
difference calculations which can be carried out to any 
desired accuracy. 

The calculation can be considerably simplified if 
only the region close to the nose is considered. The 
first few terms of the series given by Eq. (1)—for ex- 
ample, two—are then sufficient to describe the flow 
field; then only the parameters ky and k, need to be 
determined. 

As an example, the flow about the nose of the RM-10 
body, described by the equation 


vy = 0.1333 (1 — x)x 


has been analyzed at several Mach Numbers. Due to 
their negligible contribution, entropy terms have been 
neglected in this numerical work. To compare the 
solutions with data available in the literature, the rela- 
tive ratio of shock to body curvature have been plotted* 
in Fig. 3 against the combined supersonic-hypersonic 
similarity parameter in the sense of reference 5. Com- 
parison with results given in reference 5 shows that the 
agreement is indeed very good. 


* A material error in such a plot for J = 3.18 exists in reference 
1 as it was courteously pointed out to the authors by Dr. M. D. 
Van Dyke. 
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LIMITATIONS OF THE METHOD 


Clearly two fundamental assumptions are at the 
basis of the validity of any linearized analysis: (1) 
the linearized perturbation parameters must be small 
so that their second or higher order terms can be 
neglected, and (2) the derivatives of the disturbances 
must be of the same order as the disturbances them- 
selves. Failure to satisfy either condition invalidates 
the results. Unfortunately only the validity of the 
first assumption can be checked a priori. Little effort 
has been devoted to deriving a dependable way of sub- 
stantiating the actual compliance with the second as- 
sumption. However, if the perturbation derivatives 
are not finite, the flow pattern will exhibit a tendency 
to produce envelopes of characteristic lines; this fact 
might be used to set forth a criterion to determine the 
limitations of the method with respect to the second as- 
sumption. 

The case of a body shape described by the first two 
terms of the series given in Eq. (1) shall be dealt with 
first. The characteristic lines of the basic flow plus 
the disturbance flow are defined by 


(dy dx); = tan (u + ¢) = 
tan (@#+o+m+¢) (195) 


where, to the linear approximation, 


= Rix) = Rk, fil) 


uw = a+ (y — 1) 2] X (16) 


Referring to Fig. 4, if AA; and BB, are two elements 
of characteristic lines of the first family, and AA, and 
BB. are two elements of characteristic lines of the 
second family, the angles 7; and 72 between the radius 
and the characteristic lines are expressed by 


+h) 
(17) 


R ki( fo — fi) 


Envelope formation can only occur if the angles 7; and 
and 7, increase with R. Accordingly, from Eqs. (17) 
the possibility of the envelope of characteristic lines 
of the first family exists only when, for all y, 

ky (fe + fi) > 0 (18) 
while for the second family characteristic lines this 
occurs when, for all y, 

ki (fo — fi) > 0 (19) 
It has been shown that and (W) are independent 
of the sign and magnitude of body curvature and are 
only functions of the basic flow properties. Moreover, 


when the basic flow is taken as the flow about a cone, 
the following inequalities hold for all y: 


>0, AW >0, > f(y) (20) 
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Hence, for convex bodies (k; less than zero) only en- 
velopes of second family characteristics must be con- 
sidered, while for concave bodies (k; greater than 
zero) only envelopes of first family characteristics can 
occur. The region where envelopes can appear de- 
pends on the distance R from the apex. 

Call do the variation of the radius, along the char- 
acteristic line through A, corresponding to an angular 
variation dy (see Fig. +). The rate of change of this 
quantity with is Oy = A(R, and it must be 
an increasing function of R, for all y, if envelopes are 
to be avoided. Thus, 


OA OR > 0 (21) 
forall y. Considering that 

A = R tan 
and that, to the linear approximation, 


tan 72 = tan + + 
Rk, (f2 — fi) [1 + tan®* (@ 


the inequality of Eq. (21) becomes 


(fo — fi) 
— ¢+ wp) 


[tan(@ G+ = 


If the basic flow is free of envelopes, the condition for 
which envelopes of second family characteristics will 
not form is obtained as 


R< (1 ky (fo — fi) (22) 
A similar condition for the other family characteristics 
is 

R < (1/4) sn 2 (@+ 2 (fo thi) (23) 


For values of R somewhat lower than the values given 
by the equalities in Eqs. (22) and (25) the derivatives 
of the disturbance components still have finite values 
and are small. In that region the approximation given 
by Eq. (2) is then valid. 

In applying these results to the example considered 
in the previous section, it is noticed that the subject 
body has a parabolic shape with positive curvature 
and, therefore, only second family characteristics could 
formenvelopes. The limiting radi 


R = (1/4) sin2 kilfe — fi) 


~—— FIRST FAMILY CHARACTERISTICS 
—-- SECOND FAMILY CHARACTERISTICS 


Fic. 4. Characteristics of the actual flow. 
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are given in Fig. 5 for the Mach Numbers considered. 
In the region for which R (y) < R(y) the total flow 
field is certainly free of envelopes and the validity of 
the results have only to be checked against the accuracy 
of the assumption concerning the order of magnitude 
of \;. As should have been expected, the limiting 
radius becomes smaller and smaller as the combined 
supersonic-hypersonic similarity parameter increases. 
The extension of the analysis to more general shapes 
to take into account higher order terms is rather evi- 
dent. For a body described by Eq. (1), Eqs. (15) and 
(16) are replaced by the following: 


(dy dx)y. = tan(u+ ¢) = 


n 
tan E +o+ Dut 


j=1 


= (¢,/kyw = Rk fi. (W) 


j=l j=l j=1 


[sin? + (y — 1)/2] (V, 
cos 
= R’Ri fs, 5 (W) 
Hence, 
a-etwt 
R’k,( fe, 5 — fr, + tan’ 


tan 72 = tan ( 
The inequality given by Eq. (21) must still hold, and, 
if the basic flow is free of envelopes, it becomes 


sn2(@—¢+y)=0 (24) 


A similar inequality for the second family of character- 
istics is 


1 — E G+ 1) fo, +h, + 
j=l 


(25) 


The problem of determining the range of validity of 
the solution is dealt with as before. Thus, for in- 
stance, when the body is described by a cubic (7 = 2) 
the region of flow free from envelopes of either family 
is defined by the following inequalities: 


(26) 
1 — [4RRi(fo.1 + fi, 1) + GR*Ro( fo, 2 + fr, 2)] + 
sn (27) 


It is interesting to point out that envelopes of the first 
family characteristics can now appear even when the 
body has positive curvature. Indeed Eq. (26) has a 
positive root R = R* if ko(fo,2 — fi,2) > O and the 


AERONAUTICAL 


SCIENCES DECEMBER, 1957 


7 ~~ BASIC” CONE 7-7-7 


| 
Ze 


| } | 
6° 10° 14° 18° 


Fic. 5. Limiting radii for the RM-10 body with a 7°30’ angle 


basic cone at several Mach Numbers. 


polynomial expression given by that equation becomes 
negative for R > R*. 


CONCLUSION 


A rapid and accurate method to determine the flow 
field about a body of revolution has been presented. 
The method is based upon the linearized characteristics 
approach discussed in reference 7. Comparison of few 
examples with available data shows that rapidity has 
not been achieved at the expense of accuracy. An 
important useful feature of the method consists of the 
possibility of obtaining an upper limit of the validity 
of the assumption that the perturbation derivatives re- 
main finite and small throughout the flow field. 

The characteristics of the method make it worth ex- 
tending to more general fields without axial sym- 
metry. This is presently being done. 
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The Response of a Bisymmetric Aircraft to 


Small Combined Pitch, Yaw, 


and Roll 


Control Actions 


ROBERT A. DAVIS** 
The RAND Corporation 


SUMMARY 


The linearized equations of motion of a rolling bisymmetric 
aircraft are formulated, using body axes. It is shown that (a) 
the rate of change of angle of attack is a vector, (b) a linearized 
form of the downwash-lag terms is valid only if the downwash- 
roll angle be small, and (c) such a linearized form must use the 
vector relations of (a). Asa result, stability criteria are derived 
which are independent of the rolling of the reference axis system. 
The complete solution to the equations is presented for a stepped 
The complete 


pitch/yaw control action at constant roll rate. 
valid only for a 


solution is likewise given for a variable roll rate, 
simplified system in which the Magnus Effect is neglected. 


SYMBOLS 


A = incremental roll rate, rad. /sec. 

D = differential operator with respect to time 

Ey = determinant of the matrix u 

F = adjoint of the matrix u 

I = the unit matrix, or moment of inertia about the 
VY and Z axes, lb. ft. sec.” 

moment of inertia about the Y axis, Ib. ft. sec.? 


= 

K = modal matrix of the matrix u 

L, Jl, N. = moments about the Y, V, Z axes, ft. Ibs. 

ly = distance between fore and aft fin centers of 
pressure, ft 

m = mass, lb. sec.?/ft. 

P,QV,R = angular velocity components about the Y, Y, 
and Z axes, rad./sec. 

t = time, sec. 

U,V, WW = linear velocity components along the Y, }, and 


Z axes, ft./sec. 
system characteristic matrix 


u = 

t = downwash velocity, ft./sec. 

2. = force components along the Y, }, and Z axes, 
Ibs. 

y = column matrix of variables 

a,B = angles of attack in the Y-Z and Y-J planes, rad. 

a,, @,, @ = vector components of the rate of change of angle 


of attack about the Y, V, and Z axes, rad. /sec. 
= flight path angle, rad. 


Presented at the Stability and Control Problems During Large 
Maneuvers Session, Twenty-Fifth Annual Meeting, IAS, New 
York, January 28-31, 1957. 

* This paper is an abridgment of a Doctoral Thesis submitted 
to the Graduate Division of the College of Engineering in partial 
fulfillment of the requirements for the degree of Doctor of Engi- 
neering Science of New York University, October, 1955. 

** Engineer. 

7+ The author’s attention has been called to NACA TN 3737, 
“The Motions of Rolling Symmetrical Missiles Referred to a 
Body-Axis System,’’ by Robert L. Nelson, November, 1956 
This note is a valuable extension of reference 1 and should be 
referred to for amplification of both reference 1 and the present 
paper, provided, however, that the roll rate limitations and their 


cause as cited herein are kept in mind. 


A = column matrix of system inputs 

b,, 6,, 8, = control surface deflections in roll, pitch, and 
vaw, rad. 

€ = downwash angle, rad 

d, = sth root of the matrix uv, 1/sec. 


real part of \,, 1/sec. 


= 
Q = matrizant; complementary function matrix 
w = imaginary part of \,, 1/sec. 
Subscripts 
c = complementary function 
0 = initial condition 
P,Q, R, U = partial derivative with respect to ?, Q, R, and | 
s = steady state 
t = transient 
= Xx, y, component 
a, 8,4 = partial derivative with respect to a, 8, or control 
action 
a, 8 = partial derivative with respect to components 
of rates of change of a or 8 
Pa, PB = partial derivative with respect to product of roll 


rate and a or 8 


INTRODUCTION 


AND attempted to extend classi- 
cal ballistic treatments (see, for instance, refer- 
of spinning shells to bisymmetric spinning 
aircraft. Their reference axis systems differed; Nico- 
laides’ pitched and yawed but did not roll, Bolz’ did all 
Their stability criteria differed as to the role 
This paper re- 


ences 4 


three. 
played by the downwash-lag terms. 
conciles that difference by using body axes, as did Bolz, 
to derive the approximate criteria of Nicolaides in 
closed form. In the process, it is seen that a lineariza- 
tional bound on the size of the roll rate exists. Char- 
ters’ made clear the fundamental assumption under- 
lying Nicolaides’ work; for a bisymmetric configura- 
tion, the aerodynamic forces and moments are invariant 
with static roll angle. He assumed that, except for 
Magnus effects, this be true as well with roll rate. As 
shown here, this is true only if the roll rate is small. 


FORMULATION OF THE EQUATIONS OF MOTION 


General 

The linear formulation used by Jones in reference 9 is 
adapted here. The body axis system is shown in Fig. 
1. The bisymmetric aircraft is assumed in horizontal 
straight and level flight at constant speed with constant 
mass. Gravity is neglected. The variables desired, 
lV, W, Q, and R, are disturbance velocities. 


< 
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At the outset, |’ and W will be replaced by their 
‘small-angle’ equivalents, — l’6 and Va, respectively. 
They will be used again after the following discussion 
of the rates of change of a and 8 with time. 


The Downwash-Lag Phenomenon 


The Effects of Downwash—The downwash-lag effect 
exists because the air must take a finite time, /7/U, 
to travel from the front to the tail surfaces. During 
this time the free-stream angle of attack can have 
changed. The effect is then interpreted in terms of the 
instantaneous difference in local angle of attack on the 
tail surfaces over that holding for an unchanging free- 
stream angle of attack. Furthermore, it is measured 
in a plane normal to the tail surfaces. 

Note that the free-stream angle, a, and the downwash 
angle, ¢, it causes are measured in planes normal to the 
front and tail surfaces, respectively. These planes 
remain coincident during simple pitching. They do 
not during rolling. Consider, for example, the sche- 
matic tail view of a rolling bisymmetric missile in Fig. 
2. 

At time 4, = O, the relative free-stream velocity com- 
ponents are LU’, 1%, and WW). At time 4 = /7/U, they 
are lL’, V’, and IW’, and the aircraft has rolled through an 
angle P/y LU’. We wish to compute the differences in 
downwash velocities normal to the yaw and pitch tail 


surfaces at time /;. These differences are 


w, = +(dw./OW)} Vi — Vo cos (Plr/U) — 
Wy sin /U)! 

= +(dw./0W)} Wi — Wy cos (Plp/U) + 
Vy sin U)! 


(1) 


respectively. 

The downwash-lag forces and moments are directly 
proportional to these differences. As it is, however, 
Eq. (1) must be linearized. The downwash_roll- 


NOTE. 

POSITIVE FORCE,MOMENT, LINEAR AND ANGULAR VELOCITY 
COMPONENTS ARE CONSIDERED POSITIVE VECTORS ALONG 
AXES.CURVED ARROWS SHOW SENSES OF ROTARY ACTIONS 

INCLUDING THOSE OF CONTROLS. 


V/A 
A/eA/ 2 
* $/8 SERIE 
Rin] 8 Wy 
Fic. 1. Axis system and sign convention. 
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V, 
° Y 
| Wo 
Z 
t=O t=£,/U 
Fic. 2. A rolling missile. 


angle, P/; U, must be kept small. This yields 


Aw, & (dw./9W){ — WiP}r/U 
Aw, & (dw, + U 


(2) 

Clearly, then, the downwash-lag terms are not only 
proportional to the scalar derivative of the normal free- 
stream velocities but must include cross-coupled veloc- 
ity terms as well. Put another way, the classical use 
of a or 6 to describe the effect is insufficient where roll- 
ing is considered. This is emphasized by dividing Eq. 
(2) by U and getting 


—(de da)\B + Pallr/U 
+ (de da)\a PB\ lr U 


Aes 


(3) 


The Rate of Change of .\ngle of .A\ttack 
cance of the added terms in Eqs. (2) and (3) goes deeper. 
Consider them from the kinematic point of view. Rela- 
tive to inertial space at any instant, there are two angu- 
lar rates of interest. One is the space rate of body 
rotation given by 


Q = Pi + Qj + Rk (4) 


The signifi- 


where 7, j, and & are unit vectors along the x, y, and z 
axes. The other is the space rate of velocity rotation 
given by 


r= V-'Xe (5) 
where 
V = Ui+ Vi + Wk (6) 
a = a;i + a,j + a,k (7) 
a,= VR+ WOQ) 
and a, = V — WP + UR} (8) 


a, = W-— UQ+ VP 
Their difference, 
A=o-fT (9) 


is the vector rate of change of angle of attack. Its 
components are 


a, = P—(1/|V|){ VW WV-(W? + VP + 
U(WR + VQ) 

WU — (U2? + +1 9) 
V(UP + WR)} | 
a,= 
+ 


a, = Q- (1 


These are the complete expressions for the rates of 
change of angle of attack about each axis. 


They may 
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be simplified by taking into account the relative size of 
and its constancy. This yields, to the first order, 


a, = (1/U){W + VP} (11) 
a, = (-1/U){ V — WP} 


Note that these are purely kinematic relations and 
have nothing to do with either the configuration or the 
aerodynamics involved. 

Downwash Lag in the Rolling Missile—Taken to- 
gether, Eqs. (3) and (11) show that use of the scalar 
rate of change of angle of attack alone is inadequate 
from both the aerodynamic and the kinematic point of 
view. But it has not been shown that these two views 
are mutually compatible. In other words, can P be 
set equal to zero in Eq. (2) if nonrolling axes are as- 
sumed ? 

The answer is yes. The reason is that the same com- 
parison principle was used to derive Eq. (2) as was 
used by Euler! in deriving Eq. (8). Velocity differ- 
ences are measured after a small rotation of the system 
is allowed. Providing the downwash roll-angle, P/,/ V, 
is kept small, (say, 0.1 or less) 1’ and IV in nonrolling 
coordinates are identical to 1 — WP and W + VP in 
rolling, as body axes, coordinates. Thus, we have the 


inherent restriction that 


P< 01 Ir (12) 


The Equations of Motion 
The necessary alterations to the classical downwash- 
lag terms discussed above are the chief modifications 
to the equations of motion as given by Jones*. For 
convenience, the equations are written now in matrix 
form: 
(ID + uy = A(t) (13) 
P—rP = —Pyr (14) 
Eq. (13) is the system of linear first-order differential 
equations for V, I’, Q, and Rk. Eq. (14) is the adjoint 
equation for P. 
The meaning of the undefined symbols in Eqs. (13) 
and (14) is as follows: 


a —bP 0 U 
bP a —U 0 (15) 
= 
cP —d —fP 
d cP fP e 
(oee—du)<0 
UNSTABLE 7 UNSTABLE 
STABLE 
| 
| 
| | 
4 
b-f 


Fic. 3. Stability boundary for a stable nonrolling missile. 


In Eqs. (15) and (16), 


a= —(Z,/mU), e = -((M, + M,) 

b= 1+ (Zpg/mU), f 1 — 

c = (Mp,/IU) + g=Z,m 
(ZpgM,/mU?I), 

d = (M,/IU) + h 


(M; + (Z;M, mUT) 
(20) 
STABILITY CRITERIA 


General 


Kelley and McShane* and later, Murphy,’ included 
Magnus effects in their extension of the basic work of 
references 5 and 6 on spinning shells. These require no 
consideration of downwash lag. Bolz* attempted to 
account for downwash lag using body axes in extending 
the analysis to bisymmetric missiles but did not include 
the required extra terms of Eq. (2). Nicolaides cir- 
cumvented the difficulty by using nonrolling axes but 
derived only a first-order expression for the stability 
criteria. Bolz’ criteria do not agree with Nicolaides’ 
in the role assigned to downwash lag. In fact, neither 
is at once demonstrably independent of whether the 
reference axis system rolls or does not. Using body axes 
as did Bolz, the closed form equivalent of Nicolaides’ 
criteria will now be derived after the procedure of 


reference 3 and 7.* 


Derivation 
The characteristics of the free motion of the system of 


Eq. (13) are epitomized in its characteristic equation: 


Ad! + + DA+ E = 0 (21) 


where 
B=2(a + e) 
C = (a + e)? + 2(ae — dU) + (8° + f*)P? 
D = 2(a + e)(ae — dU) + 2P? + af? + (f+ 
E = (ae — dU — fbP?*)? + (be + af + cl’)*P? 
(22) 
The roots of Eq. (21) are 
= ®) + + 9) 
Ae = (E — R) + Uy — 9)! (23 
As = (E+ — i(n + 9) 
Ay = (E — R) — Un — F) 


* An unidentified referee has informed the author that ‘‘no 
attempt was made in reference 2 to take into account the wing- 
tail downwash effects for a rolling missile.’"’ This was certainly 
not clear from the text of reference 2. 


= 


907 
~ PA /u A(t) = —gé,, g6,, (17) 
P, = —LZ,/Lp (1S) 
j 
Is 
ee 
(6) 
(7) 
4 
= 
= 


908 JOURNAL OF THE AERONAUTICAL 
pé 
4 
| (ae—dU)>0 
! 
1| STABLE { 
if | 
3 | 
| a 
UNSTABLE UNSTABLE 
| | 
| | 
| | 
4 
‘b-f 


Fic. 4+. Stability boundary for an unstable nonrolling missile. 


where &= —(1/2)(a + e) (24) 
n= +(12)(64+/f/)P (25) 
= [—w? — — + 
P2[(1/2)(a — e)(b — f) — + 
[—w? — (1 4)(b — (26) 
— (1 406 — + 
P2[(1/2)(a — e)(b — f) — 
[—aw? — (1 4)(b — (27) 
wy? = ae — dl — (1/4)(a + e)? (28) 


Analysis 
From Eqs. (23) the free motion will be dynamically 
stable if and only if 


(29) 


In order to study this further, it is necessary first to 
consider the stability of the nonrolling aircraft. 

Eq. (2S) is the damped natural frequency of the non- 
rolling aircraft provided that 


(ae — dl’) > (1 4)(a + e)? (30) 


This can be seen by setting P = Oin Eqs. (26) and (27). 
Thus, if Eq. (30) holds, then 


0: >0 (31) 


If Eq. (30) does not hold, thén 


w 
to 


—w?>0; ¢ =0 (: 
It follows that if 

(ae —dU)> (1 4)(a +e)? > O (33) 
or (1/4)(a + e)? > (ae — dU) > 0 (34) 


the free motion of the nonrolling missile is dynamically 
stable and either oscillatory, Eq. (33), or nonoscilla- 
tory, Eq. (34). If 


(ae —dl) <0 (35) 


it is dynamically unstable. 
Returning to Eq. (29), we obtain by substitution 


Q 
Wes = QusVo + § 
= 0315 Vo 
= Vo ¢ 
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from Eqs. (24) and (26), 


Pb — (6b — f)P (a — e) X 
(b ae} < (a + e)*(ae — dU) (36 
This is the basic stability criterion for bisymmetric 
If it is satisfied, the free 
If not, then it is dy. 


missiles at constant roll rate. 
motion is dynamically stable. 
namically unstable. 

Consider the case in which the nonrolling aircraft js 


dynamically stable—that is, when 


(ae —dl’)> 0 (37 


This case is shown in Fig. 3. Here, the roll rate at 
which Eq. (36) is not satisfied is plotted as a function 


of the parameter cl’ (6 — /). For the range 
—e< (38 


the motion is dynamically stable for any value of P* 
Outside this range, increasing ?° degrades stability 
until the free motion is unstable. 

When the nonrolling bisymmetric aircraft is dynami- 


cally unstable—that is, when Eq. (35) holds—it re- 


mains so when rolled except within the range of Eq. | 


(38). 
possible to stabilize the missile. 

The criterion of Eq. (36) is derivable in terms of the 
such a 


If rolled sufficiently fast there, it is theoretically 
This is shown in Fig. 4, 


coordinates of a nonrolling axis system. In 
system, the terms } and f of Eq. (36) as defined in Eq. 
(20) would each be reduced by unity. Their sum 
in Eq. (36) is reduced by P which is as it should be. 

Finally, it is seen that if the Magnus terms / and ¢ 
are neglected, Eq. (36) is satisfied identically. This 
implies that within the limits of this analysis only the 
Magnus effects cause dynamic instability of a rolling 
bisymmetric aircraft. 


STABLE SOLUTIONS AT CONSTANT ROLL RATE 
The complete solution to Eq. (13) is, in matrix form," 


y(t) = “yy + [I — JA (39 


The Complementary Function 

For comparative brevity, half of the complete ex- 
pansion of the first term in Eq. (39) is given below. 
The subscript ( ), appearing there must take the values 
of one and two alternately according to the following 
scheme. 
(25) can be written as 


The general expression for the roots in Eqs. 


= Ms (40) 
where ifp=1;q¢q = 3, thens = 1 
and ifp=2;q =4, thens = 2 


The expressions for the free motion of each variable 
may be written as 


us Vo — 01, Wo + 
+ 


(41) 
243, Ro 


31s Vo + 243,00 + 233,20 
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| 
| 
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ae 
} 
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dl (36 


symmetric 
the free 
it is dy. 


aircraft js 


(94 
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function 


(38 
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Qn, = + = cos (wf + Yu; Qa, = —Qn, tan (wf + Yn 
= + cos + Qa, = tan (wd + va 1” 
= R2(F13 cos (at + Yis,); = —Qis, tan (wl + 
Qs3, = + cos (wt + Was); Osx, = tan (wt + Yss 


In Eqs. (41) Vo, Wo, Qo, and Ry are the initial conditions. The remainder of the terms in Eqs. (42) are defined as 


ws: 


R(Smns) = R(Frns)R(ds) + (ds) (43 
(Finns) = — (ds) 
where 

= P\—pblb + + ble — a)| — blaf + be + 

= + ble — a)) — + f)P + (b+ + bfP?) + e[eU + ble 

R( Fix.) = UP\ulb +f) + (cU + be + 

= UPw,(b + f) 

R(Fy) = + — e(b + fieP + d(cU + be + af)t $4) 

= + ficP + + + cP(cU + be + 

R( Fy.) = + — + fla — e)] — fP(cU + be + af); 

I( = + fla — — + f/fP + (b+ f)\(dU + fbP?) + + fla — 


and K and J signify the real and imaginary parts, respectively, of F,,,., 


Furthermore, in Eqs. (45), 


d, = 2(Fys + F335) (45) 


Finally, in Eqs. (42), 


The Particular Integral 
For convenience, the expansion of the second term of Eq. (39) is handled in t vo parts. 
The Steady-State Terms—-These are given by the expansion of the term w~'A: 
= (] { )6, + ( Fyoog + 
Wo = E»)} + (Frog + Fy 39h 6, + (Frog (47 
7) 


QO, = (1 —(Faog — Fowh)6, + (Frog + 
— (1 + ( + 6, Fyyoh 
Here, 
= Fray (ae — dU)e + (af + 
Frio = —Foa = (be + cU)eP + (dU + fbP?)fP 
Fi39 = Foy = (be + af + cU)UP 
Figo = — Foo = (—ae + dU 4+ OfP?)U (48) 
Fy) = Fy = | — (be a af)d + (fbP? — ae)c|P 
Fy = — (ae — dal’ — bfP*)d (af + be + cl)cP? 
F320 = = (ae — dl ja + eb)bP? 
= — = [(af cl)a + (dU + bfP?)b|P 
and ky = (ae — dU — jor* + (be + af + cl’)*P? (49) 


The Transient Terms—These are essentially a composite of the complementary and the steady-state terms, as is 
evident from Eq. (39). Thus, we define the matrix 
T = (1/ Eo) (50) 
where the elements of the matrix F, are given by Eqs. (48). The elements of the matrix 2 are given by Eqs. (42). 
Under the same subscript proviso holding for the writing of the complementary function then, the transient terms 


may be written as 


- = (Trig — (The. + Tre 
0, = — — (Tag + 


R,, = + T33,/2)6, (Ta1.2 


4 
‘ 
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CONTROL INPUT 


In this section the Magnus effects are neglected. This 
is necessary to derive the desired matrix solution. It, of 
course, ensures dynamic stability as is clear from Eq. 
(36). Beyond the mathematical convenience, how- 
ever, neglect of the Magnus effects on a rolling missile 
should be effected with great caution. This has been 
amply emphasized by Nicholaides.! One safeguard is 
that we assume the aircraft to be rolling sufficiently 
slowly and to experience the effects of a small stepped 
roll control input. 

Proceeding, therefore, the Magnus terms in c and } 
are set equal to zero in Eq. (15). Furthermore, / is 
set equal to unity under the assumption that /, / is 
negligible. The roll rate is now a function of time. 
The roll inputs are assumed to be introduced simul- 
taneously with the pitch and yaw control inputs. 

Eq. (15) becomes 


a — P(t) 0 l 

u(t) P(t) a —l 0 (52) 
—d — P(t) 
d 0 PO) e 


The roll rate is given by the solution to Eq. (14). 
(53) 
A=P,-—P, (54) 


where 


Substitution of Eq. (53) in Eq. (52) permits Eq. (13) 
to be rewritten as 


[ID + u, + IAe™|y = A (55) 

where 0 0 
] 0 0 0 
0 0 -1 

0 0 


and where “, is now a matrix of constants containing 
P, instead of P(t) asin Eq. (52). Eq. (55) is the system 
of equations to solve. 

The solution of Eq. (55) is most conveniently ex- 
pressed in terms of its matrizant operator. Thus, 


y(t) = Q(—Uu, + 


Q(—u, — Q-(—u, — IL.ledt-A (57) 
0 
where 
Q(—u, — = Q(—u,)2(—T Ae”) 
= e “(I cos (A /r)(e — 1) 
II sin (.1 /r)(e’"’ — 1)] (58) 
It can be shown that the matrices II and e~“" are 


commutable. Furthermore, the matrix e~““ will be 
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recognized by its components already given for nop. 


zero b and c in Eqs. (42). 


The complementary solution can be written in closed f 


form immediately as 


ve = [I cos /r)(e" — 1) 
II sin (1 — 1) y(to) (59 


The particular integral must be expanded in series 


however. These terms are given by 


A fu, ! 
Limo r ’rK=0 
| "4 n 1 
(- II ) + IK) (60 


n=0 Y K=0 


The first term of Eq. (60) corresponds to the “‘steady- 
state’ terms. It cannot be constant as with a constant 
roll rate, however, but must change as the roll transient 
proceeds. As can be seen, it eventually approaches 
that value holding for a constant roll rate, P,. 

The second term of Eq. (60) represents the transient 
This is seen to be the product of a constant matrix 
and the transient factor written for a constant rol] 
rate, P.. 

The crux of the above solution lies within the com- 
mutability property of the matrices e-““ and II. 
permits the development of the matrizant in the closed 
form of Eq. (58). Such cannot be done unless the 
Magnus terms in 6} and c are set equal to zero and f is 
set equal to unity, as examination of Eq. (15) shows. 
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The Effect of Molecular Vibration on Recovery 
Temperature in Plane Couette Flow; 


Hsun-Tiao Yang* 


Institute for Fluid Dynamics and Applied Mathematics, 
University of Maryland, College Park, Maryland 


May 9, 1957 


O' A POLYATOMIC GAS, the translational as well as rotational 
degree of freedom is always active, while the inert degrees 
of freedom—namely, vibration, dissociation, and ionization—are 
excited only at relatively high temperatures. For example, for 
diatomic gases such as oxygen, nitrogen, or air, these inert de- 
grees of freedom become important at temperatures of about 
600°K., 3,000°K., and 10,000°K., respectively. The purpose 
of the present note is to study the effect of molecular vibration 
on the recovery temperature in plane Couette flow with the 
temperature range of interest between 600°K. and 3,000°K. 
The effects of dissociation and ionization have recently been 
studied by Liepmann and Bleviss,': ? who, however, neglected the 
effect of molecular vibration. 

The plane Couette flow consists of the fluid flow between two 
infinite, parallel plates at a distance 6 apart with the lower plate 
stationary at a temperature 7), and the upper plate moving with 
a velocity U in its own plane at a temperature 7..** It is as- 
sumed that sufficient time (larger than relaxation time) has 
elapsed since the start of motion of the upper plate so that the 
inert degrees of freedom are already in equilibrium; therefore the 
only independent variable is the ordinate y normal to the plates. 


+ This work was carried out under Contract AF 18(600)993, sponsored 
by the Office of Scientific Research, Air Research and Development Com 
mand, and under the supervision of Professor Shih-I Pai. 

* Research Associate 

** In order to compare with boundary-layer flow, the notations of Liep 
mann and Bleviss are followed—i.e., subscript w refers to the lower plate 
and subscript © refers to the upper plate 
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Fic. 1. Effect of molecular vibration on recovery temperature 
in plane Couette flow. 


The continuity equation and the boundary conditions on both 
plates require the normal component of flow velocity to vanish— 
i.e.,v = 0. The momentum and energy equations become 


dr/dy = 0 (1) 
(d/dy)(p +0) =0 (2) 
(d/dy) (ru — q) =0 (3) 
where 
p = pressure of the gas 
uw = flow velocity in the x direction 
7 = shearing stress 
o = normal stress in the y direction 
g = heat flux in the y direction 


It is further assumed that the flow may be considered as con- 
tinuum, so that the following Navier-Stokes approximation are 


valid 
rT = w(du dy), (4) 
= —NdT/dy) = —(1/Pr)yc,(dT/dy) 
where 
7 = temperature of the gas 
aw = coefficient of viscosity 
A = coefficient of heat conduction 
Pr = Prandtl Number 
¢,» = specific heat at constant pressure 


= R + (3/2)R + (2/2)R + 
1(07/27)?/[sinh (@7/27)|2{ R  (5)* 


im which the second term is the contribution to internal energy 
due to translation, the third due to rotation, and the fourth due 
to vibration, R being the gas constant, and 


0, = hw/2rk (6) 


where i; = Planck’s constant, k = Boltzmann’s constant, and 
w = angular frequeney of molecular vibration, 
From Eqs. (2) and (4) it is seen that 


p = const. (7) 

and the enthalpy, /, is related to temperature by 
dh = ¢,dT (8? 
Since the Prandtl Number is nearly constant for the present 


* See for instance, Eq. (14.41) p. 369. reference 2. 
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temperature range, Eq. (3) may be integrated twice as in the work 
of Liepmann and Bleviss using Eqs. (4) and (8) to give 
(u?/2) + [(h — = —(qw/tw)u 
where 7, is the shearing stress at the lower plate, and y., is th, 
heat flux from the lower plate to gas. 
With Eqs. (5) and (8), Eq. (9) vields 


(u?/2) + (7/2) (1/Pr) (T — Tw) + 
(R6,/2) [coth (0,./2T) — coth (6,/2Tw)] = —(qu/tulu (1 


The most interesting quantity is the recovery temperature, 7) 
which is defined as temperature at the lower plate being hea 
insulated—i.e., 7y = Tw for g, = 0. From Eq. (10) we obtain 
with « = U, T = T.~ at the upper plate. 


Tr + (0,/7) coth (6,/2Tr) = Tx + (0,/7) coth (@,/2T 2) + 
ToPr(0.2)Mo? (1 


where the last term may be so written, provided the upper wal 
iscool (Ts < 600°K.). It is knownt for Os that 6, = 2230°K. an 
for No, 0, = 8340°K. The recovery temperatures with molecular 
vibration taken into account are obtained for oxygen and nj 
trogen at Pr = 0.70, To = 300°K., as shown in Fig. 1. B; 
neglecting molecular vibration, one has the usual recovery tem- 
perature as, with y = 1-4, 


Tr = + Pr{(y — 1)/2 (12 


which is also plotted in Fig. 1. It is seen that the decrease of re. 
covery temperature due to absorption of internal energy by the 
molecular vibration is rather appreciable. 
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Heat Transfer in a Laminar Boundary Layer 
From a Surface Having a Temperature 
Distribution: 


B. N. Pridmore Brown** 

Research Assistant, Institute of Aerophysics, University of Toronto 
Toronto, Canada 

May 17, 1957 


SYMBOLS 


= constant derived from Sutherland's formula 

kj = free-stream thermal conductivity for air 

Po = free-stream stagnation pressure 

Pe = pressure measured by boundary-layer pitot probe 
q = rate of heat transfer 

Re = Reynolds Number 

Ti = free-stream static temperature 

Ty = adiabatic recovery temperature 

Tu = wall temperature 

t*(x) = [Tw(x) — Te)/Ti 

t*(xa) = points at which Tchebichef polynomial is fitted 
x = distance from model leading edge 


DISCUSSION 
H™ TRANSFER and pressure measurements were made in 
a laminar boundary layer generated on the outer surface 
of an unchoked hollow cylinder aligned axially with a Mach 
2.5 flow. Nilvar wire ribbon elements wound onto the cylinder 
provided a means of measuring and controlling average surface 


Tt See for instance, p. 405, reference 2. 
t This is an abridged account of the investigation reported in reference 1 
The work was supported by a grant from the Defence Research Board o! 


Canada 
** Now at Convair, Fort Worth, Texas 
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temperatures over 1/2-in. intervals. Each 1/2-in. element, con- 
sisting of a number of turns, was connected to a Wheatstone 
bridge as one arm. The remaining resistances in each bridge 
were of enough capacity that they could carry the heating current 
for the Nilvar arm. Power input to heating elements could be 
found from bridge input voltages and the resistances of the ele- 
ments. The resistance of a Nilvar element, and hence its tem- 
perature could be calculated from bridge input and unbalance 
voltages together with the thermal coefficient of resistivity for 
Nilvar. The model is shown in Fig. 1. 

Transition from a laminar to a turbulent boundary layer, ob- 
served by schlieren and confirmed by pressure traverses, occurred 
at about Re = 1.5 X 106 on the unheated model. With the sur- 
face temperature distribution of Fig. 2, transition moved forward 
to Re = 0.8 X 10°. 

Measurements were made with several surface temperature 
distributions, of which Fig. 2 is an example. For comparison 
with theory the temperature distribution was fitted by fourth 
and sixth degree Tchebichef polynomials. Fig. 3 shows experi- 
mental heat-transfer results compared with the general theory? 
and with the theory as specialized to an isothermal surface. 
There is qualitative agreement with theory, but the experimental 
results are about 80 per cent high, in agreement with those of 
Slack’ at the forward part of his model. 

Fig. 4 shows a boundary-layer pressure profile at 2.3 in. from 
the leading edge and with the temperature distribution of Fig. 
2. The experimental curve for an unheated surface is plotted for 
comparison. Abscissas of the experimental points on the un- 
heated model were about 8.5 per cent higher than those of the 
theoretical curve,? probably due to a probe effect.‘ The ab- 
scissas for the theoretical curve on the heated model are in- 
creased by the same percentage. The experimental pressure 
profile on the heated model shows a characteristically laminar 
shape, but the boundary layer is thicker than theory predicts. 
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The Aerodynamic Force Coefficients of Yawed 
Slender Configurations at High Mach 
Numbers 


Leon Trilling and James W. Clark* 

Assistant Professor of Aeronautical Engineering, Massachusetts 
Institute of Technology, Cambridge, Mass., and Research 
Engineer, United Aircraft Corporation, East Hartford, Conn., 
Respectively 

August 2, 1957 


HIS NOTE PRESENTS a simple method of evaluating the aero- 
dynamic force coefficients on a highly yawed slender body 
at high Mach Numbers. The basic idea is that if the Mach 
Number is of order 5 to 7 and the angle of yaw is of order 8° to 
12°, the effective cross-flow Mach Number is of order 0.75 to 


* 


* This note summarizes some results of a thesis submitted by Mr. J. W. 
Clark in partial fulfillment of the requirements for obtaining a Master's 
degree in Aeronautical Engineering. The authors wish to thank Mr. W. 
Taylor of United Aircraft Corporation for supplying the models on which 
the present measurements were taken. 
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1.5. The splitting of the total flow field into an axial flow and a 
crossflow is still legitimate, but the analysis of the cross flow is 
difficult since it is the plane transonic flow past a blunt body with 
turbulent separation on the leeward side. 

Since no simple technique now available permits the calcula- 
tion of such a flow, the present investigation consisted of meas- 
uring the pressure distribution on the windward side of appro- 
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priate blunt bodies in the transonic wind tunnel to provide the 
data needed for the evaluation of the side force coefficient in the 
hypersonic regime. 

The Massachusetts Institute of Technology facility where the 
tests were conducted is of the intermittent blowdown type wit; 
slotted walls; the test section is octagonal, 22 in. across the 
flats. Fig. 1 shows the basic steel structure of the model and 
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Variation of L/D with a for the NACA RM-10 body at 
a Mach Number of 6.9. 
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the location of the pressure taps. The model was mounted 
vertically from the test section floor and terminated with an end 
plate; subsequent tests verified that the flow was essentially two 
dimensiozial at the taps. In order to achieve other values of the 
body diameter to wing span ratio (D/d), shells made of laminated 
birch, each having its own pressure orifices and tubing, were 
screwed into place over the basic structure. Thus, the model 
frontal area remained constant while the ratio D/b increased 
from 0.5 to 1.0 (circular cylinder). The Mach Number varied 
from about 0.68 to 1.15, with corresponding Reynolds Numbers 
of 800,000 to 1,150,000 (based on the 2-in. width). 

Fig. 2 presents the more important data in graphical form. 
It will be noticed that the drag rise is gradual. Variable Mach- 
Number runs were made, and the continuously recorded pres- 
sure data showed no peaks or discontinuities. 

Several applications of this method have been made to various 
missile configurations, including the NACA RM-10. In general, 
it has been found that for bodies with wings or fins the skin fric- 
tion drag is small compared with the magnitude of the pressure 
forces involved, and neglecting it has little effect on the resultant 
L/D curves. However, for bodies of revolution, skin friction 
must be taken into account. 

Fig. 3 shows comparison of the present method with experi- 
mental data on the RM-10 body (a parabolic body of revolution 
with no fins!) at Mach 6.9. For this case, a skin friction drag 
coefficient of 0.03 (based on maximum cross-sectional area) was 
used. Supersonic cross-flow coefficients were estimated from 
the data and the trends indicated by reference 2. 

For the configurations treated in the investigation, this semi- 
empirical method of approach has proved very satisfactory. 
Indications are that the procedure may be extended to other 
shapes, such as a cruciform missile, with equal success. Similar 
coefficients may be obtained for a missile in roll attitude alsc. 
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Note on Turbulent Friction Factor of Flow 
Through Narrow Annuli 


L. N. Tao 

Assistant Professor, Department of Mechanics, Illinois Institute of 
Technology, Chicago, III 

August 2, 1957 


N MANY ENGINEERING APPLICATIONS there is a great interest 
I centered on the problem of flow through an annulus of fine 
clearance. An analytical approach parallel to that of the pipe 
flow is given in this note and the result shows an excellent agree- 
ment with the experiments. 

In laminar region the known solution is given, f = 96/Re with 


f = (4A/S) (AP/p) (2/V?) = (d/l) (AP/p) (2/V?) 
and 
Re = dVp/u 


where f is the friction factor, A the cross-sectional area, S the 
wetted surface, AP the pressure difference across the annular 
channel of length /, p the density, d the diametrial clearance, V 
the mean velocity in the annulus, and Re the Reynolds Number 
with the hydraulic radius as the linear dimension. 

In the turbulent flow Rothfus! and Knudsen-Katz* have found 
that the maximum velocity attainable in an annulus is at the same 
location as that of the laminar flow—i.e., 
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where 7; and r2 are the inner and outer radii of the annulus, re- 
spectively, and c the radial clearance. In the case of an annulus 
of fine clearance ¢/r; is a very small quantity and the above equa- 
tion may be written as 


= 1 + (1/2) (c/n) 


This means that the maximum velocity exists at the mid-section 
of the channel. 
Based on the Blasius’ law of friction 


To = 0.0228 
and the one-seventh power velocity profile in the turbulent flow 
u = u,(y/s)!/7 


where 7, is the shearing stress at the wall, 6 the boundary-layer 
thickness, and u, the velocity outside the boundary layer, we 
obtain the mean velocity within the boundary layer 


V = O.875u, 
and 
= 0.0288p V2(u/ Vip)®-* 


Now we let the boundary-layer thickness be equal to ¢/2 similar 
to a fully turbulent flow in a pipe, this corresponds to the concept 
that two boundary layers touch each other at the mid-section of 
the annulus. Equating the pressure force and shearing force of 
a steady state, we obtain f = 0.326(Re)~*-®. This result is 
compared with the Nootbaar-Kintner’s experiment 
plotted in Fig. 1 and do show an excellent agreement. For higher 
Reynolds Number the present note can be extended by using a 
different velocity profile. 
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On the Reiner-Taylor-Saffman Dilemma 


M. Z. v. Krzywoblocki 
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_ HAS OBSERVED a high pressure amounting to half an 
atmosphere at the center of the airspace between two 
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parallel plates when they are very close together and one of them 
rotates. Reiner interprets this experimental result as being 
due to a non-Newtonian or viscoelastic property of air—i.e., 
This 
manifests itself by the presence of normal stresses in the fluid 
in the direction of the velocity. Taylor and Saffman attempt 
to show that the phenomenon is due to the geometrical and 
dynamical imperfections in the experimental apparatus (in par- 
ticular, two cases are considered: a small error in the perpen- 
dicularity of either of the plates to the axis of rotation and a small 
vibration of the rotor parallel to its axis) and to the compressi- 
bility of the air (when the fluid is incompressible, the mean 
value of the excess pressure is found to be zero). 

The remarks below are of two kinds. First refer to some par- 
ticular items used by Taylor-Saffman. In trying to disprove 
Reiner’s hypothesis (non-Newtonian properties may be described 
by the second approximation to the relationship between the 
stress and rate of strain tensors) they use the results of the 
kinetic theory of gases in the Chapman and Cowling form. As 
is well-known, nothiag is more doubtful than the validity of that 
derivation, where the higher-order approximations, which may 
contribute never calculated. 
Moreover, the very deep and sensitive fundamentals of the 
kinetic theory of gases were never clarified. 

A lot was written on that subject without any substantial re- 
sult, and consequently, it seems that this item in the Taylor- 
Saffman discussion cannot be accepted seriously. 

The second kind of remark refers to the fundamentals of our 
approach to the description of phenomena of nature. The explana- 
tion offered by Taylor-Saffman based on the compressibility and 
analogy with the well-known action of lubricated bearings seems 
to be reasonable and acceptable. The question whether all the 
details of the phenomenon can be explained on tie basis of the 
They wrote, 


fluid has an elastic resistance to shear as well as viscosity. 


some lower-order terms, were 


Navier-Stokes equations only is an open one 
“. . . in the study of the mechanics of air-lubricated bearings 
compressibility may have importance in cases where fluid is 
dragged into a narrow space... . Cannot normal stresses be 
responsible for this dragging of air? They are responsible for 
the dragging of water up to the highest points of the highest 
trees. 

In our approach to the description of nature’s phenomena we 


” 


may distinguish three regions: 

(a) macroscopic phenomena in macroscopic domains tested 
by macroscopic instrumentation, 

(b) macroscopic phenomena in microscopic domains (or at 
least having microscopic dimensions in one direction as in the 
present case) tested by macroscopic instrumentation (or by sub- 
microscopic instrumentation ), 

(c) submicroscopic phenomena in submicroscopic domains (on 
the scale of nuclear phenomena) tested by submicroscopic instru- 
mentation. 

In each of these domains we try to describe the phenomena by 
means of some analysis—equations. Since we do not know the 
characteristics properties of the medium, we usually put forward 
some hypotheses on the nature of these properties. But the re- 
gion of validity of these hypotheses absolutely cannot be extrap- 
olated outside the region in which they were verified experi- 
mentally. This was the fact with the Newton Law breaking 
down in (c) and the origin of quantum mechanics. The Navier- 
Stokes equations were primarily proposed for and verified in (a) 
but up to now nobody verified them in (b) (sublayer in the 
boundary layer, say) since nobody possesses instrumentation 
adequate enough to test any hypotheses in (b). As a matter of 
fact we may be never able to verify experimentally phenomena 
in (b), since they and the region (b) may be partly affected and 
overlapped by (c). 

In this case the Heisenberg principle of uncertainty will be 
valid and we shall never be able to construct instrumentation 
(a) adequate for (b), unless we shall use the tool of (c). It is easy 
to show that analysis (equations) valid in (a) may break down 
in (b). Any function, or constant (like u or \) can be decom- 
posed by means of harmonic analysis into infinitely many 
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functions. 
(a) are not simple functions as given to us by the previous ce 
tury; they appear only as simple when tested in (a). But the 
may be composed of some more refined properties (function; 
unknown at the present time which in (a) do not demonstra, f 
their individualities. But in (b) due to the vicinity of other bog 
ies (discs), intermolecular and nuclear forces, and maybe othe 
reasons, some of them may demonstrate their individualities jy 
way different from their behavior in (a). 
preted as a change in value of w» in (b) or that some proper 
manifests its importance in (b) and unimportance in (a). Am 
how, the physical meaning of » both in continuum and kineti 
theory of gases is one of the most crude definitions in physi 
(particularly when modern physics enters the picture). 

Thus, the statement that the Navier-Stokes equations myy 
describe accurately and always the conditions in (b) may 
wrong in itself. It seems that Reiner, guided by intuition, fee 
that they are not quite justified in (b) and has proposed som: q 

a 


It is possible that the characteristics properties ; 


This may be inte 


thing else. It seems incorrect to state that there is one set 
equations which describes adequately (a), (b), and(c). We kno 
today that (c) can be described only by means of quantun f 
mechanics. Whether it is better to reach (b) from (a) or from (¢ 
this question cannot be answered today. I am rather of the opin- 
ion that the sure and safe way to investigate (b) is to start frop 
(c) and in the period of the last several years in many talks I have 
called attention to the necessity of deep investigation of (¢ 
(quantum hydrodynamics). The range (c) already has give: 
some fair answers in such complicated phenomena like thos 
referring to Helium I and IT. 

Fluid dynamics in its history was compelled to use as a to 
results of other fields—continuum, kinetic theory, etc. The im 
pact of modern physics upon other fields is enormous and not yet 
clearly seen. But it is very clear to me that fluid dynamics, due 
to phenomena like dissociation, ionization, and those of electro- 
magneto-hydrodynamics, must begin to use quantum-mechaniec: 
In this case (b) may be investigated on the 
way passing from (c) to (a). Theoretically, such a passage from 
(c) to (a) is possible and was already proposed. But it will be 
a long time before we shall be able to calculate all the details of 
both physical and engineering nature and importance 


results as a tool. 
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Note on ‘‘Dissociation Effects in Hypersonic 
Viscous Flows”’ 


Denis J. Zigrang 

Missile Development Division, North American Aviation, Inc 
Downey, California 

August 2, 1957 


gemma AND INTERESTING paper by Kuo! formulated equ 
tions applicable to a boundary layer in dissociative equi 
librium. During the development, the energy equation was 
written in terms of Prandtl (= c,u/dA) and Lewis (=pDe,/) 
Numbers which were later assumed constant. This note is a dis 
cussion of that Prandtl Number which differs somewhat from the 
Prandtl Number said by Kuo to be “ . . . insensitive to change of 
temperature....”’ 

The energy equation was first written as Eq. (4.3) in reference 
las 
pu(0H/dx) + px(0H/dy) = u(dp/dx) + (0/dy) [A(OT/dy) + 

(Ha — + w(Ou/dy” 


where H includes the dissociation energy; \, O, and « are the 
thermal conductivity, the binary diffusion coefficient, and the 
mass fraction of dissociated molecules, respectively; the term 
(H4 — H4,) is proportional to the difference between the sum 
of the chemical and sensible heat content of the atom and the 
sensible heat content of the molecule; and the last term in the 
bracketed portion of the energy equation, (H4 — H4,)pD(0«/09), 
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js an expression of the amount of thermal energy transferred by 
diffusion of the chemically active species along the concentration 
gradient. Thus, the thermal conductivity used in the energy 
equation does not include a portion of the dissociative equilibrium 
contribution, namely that of the diffusional mechanism. 

An isobaric heat capacity which includes the effect of dissoci- 


ition was introduced and the bracketed term of the energy equa- 


tion rewritten in Eq. (4.6)' as 
((u/Pr) {1 + ((Ha — Ha,)/c,|Le(de dT)} (O0H/dy)) 


where P7( = ¢,m/%) and Le( = pOc,/d) are, respectively, the 
Prandtl and Lewis Numbers This Prandtl] Number 
which is not consistently defined with respect to the effect of 
dissociative equilibrium was then regarded as “‘ . . . insensitive to 
.”’ based on an argument? that the 
varies very nearly 


” 


change of temperature... 

thermal conductivity coefficient A... 
as Cy, Cr being the heat capacity at constant volume . 
However, examination of Eq. (11-20)? reveals that the thermal 
conductivity coefficient which varies as ¢,u clearly includes the 
diffusional mechanism. A Prandtl Number insensitive to 
changes in temperature by the argument of reference 2 must 
therefore be written as 


Pr = IX + (H4 H4.)pOQ(de/dT )| 


It can be shown by algebraic manipulation that this redefini- 
tion of the Prandtl Number based on use of a thermal conductiv- 
ity which includes the diffusional mechanism necessitates the 
substitution of unity for the defined constant o in Eq. (5.6).! 
It was concluded! that the ratio of heat fluxes at the surface of a 
cooled flat plate for the dissociative equilibrium case and the case 
of the undissociated gas is sensibly proportional to (up/p'p®)!/? 
where »p and p refer to the undissociated gas. That con- 
clusion does not appear to be affected by this proposed change. 


REFERENCES 


Kuo, Y. H., Dissociation Effects in Hypersonic Viscous Flows, Journal of 
the Aeronautical Sciences, Vol. 24, No. 5, pp. 345-350, May, 1957 
(Ed.), High Speed Aerodynamics and Jet Propul 


2 Rossini, Frederick D 
411; Princeton 


sion, Vol. 1, Thermodynamics and Physics of Matter, p. 


University Press, 1955 


Collision of Plane Shock Waves With Wire 
Screens{ 


W. J. Franks* and J. Gordon Hall** 

Convair, a Division of General Dynamics Corporation, 
Fort Worth Division, Fort Worth, Tex., and University of Toronto, 
Institute of Aerophysics, Toronto, Canada, Respectively 


August 5, 1957 


wu EXPERIMENTAL RESULTS for collision of plane shocks 
with screens or grids have been reported,': 2 * no attempt 
at theoretical prediction of the resulting transmitted and re- 
flected shock strengths seems to have been made. The purpose 
of this note is to mention shock-screen experiments done at 
University of Toronto Institute of Aerophysics in 1954-55 and 
indicate a simple analysis used which accurately predicted the 
shock strengths measured. 

Experimental shock-tube studiest were made of collision of 
shocks in air with 30- and 8-mesh wire screens having 62 per cent 
and 50 per cent blockage, respectively. Shock pressure ratios 
used (12 maximum) provided both choked and unchoked screen 
flows. Reynolds Numbers based on screen wire diameters and 


+ The authors wish to acknowledge helpful discussions with Dr. I. I. 
Glass and financial assistance received from the Defence Research Board 
of Canada 

* Propulsion Engineer. 


** Research Associate and Assistant Professor of Aeronautical Engi- 


neering. 
} Initial experiments were done by I. J. Billington 
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Fic. 1. Idealized distance-time (x, ¢) diagram of shock-screen 
collision. State (4) and Sc present only with choked flow through 


screen. 


flow immediately upstream of screens ranged from 301 to 1,691 
for 30-mesh and from 3,040 to 8,740 for 8-mesh. Wave-speed 
drum-camera schlieren and interferometric studies (in 3 X 3 in 
and 2 X 7 in. tubes, respectively) determined the wave system 
following collision and provided sufficient data for quantitative 
evaluation of the entire flow field. 

Fig. 1 shows an idealization of the observed wave system. 
With sufficiently weak incident shock the screens were unchoked, 
in which case region (4) and shock Sc were absent. As incident 
shock pressure ratio Ps was increased, 1/, = wu;/a; increased 
sharply until eventually choking occurred at 1/; = Micr. With 
further increase of Ps, 1/7; remained essentially constant at 


be 
/e 


Wire Diam. 0.036 In. 


© £xperimenta/ Ponts 
——— Mean Expt. Curve 
Mscr = 0.306 


Based On 50% B/ock., 
x= 14 


| O-Mesh Screen 50% Blockage 
| 


fe, = Pi 


Mach Number ; immediately upstream of screen vs 


Fic. 2. 
incident shock pressure ratio Ps. 


917 

| 
| 
5 ! 
(2) 

ShocK a 

| Transmitted 

6) / Shock 

Lncident 4 
Shock ww 

Ry | 


G18 JOURNAL OF THE AERONAUTICAL SCIENCES 


8-Mesh Screen $0% Blocnage 
Wire Diam. 0.036 tn. 
© perimenta/ Points 
Theory 4 
6 
B, Be 
5 35 


/ / 
2 3 4 6 7 é 
P 


6/ 
Fic. 3. Transmitted shock pressure ratio P2; and reflected shock 
pressure ratio P55 vs. incident shock pressure ratio P¢. 


Mscr. The flow expanded through the screen from M;cr to M,> 
1 in supersonic state (4) which was terminated by shock Se facing 
upstream but drifting slowly downstream. 

In order to theoretically predict transmitted and reflected 
shock pressure ratios for given Ps and screen blockage, the wave 
system of Fig. 1 was assumed—i.e., one-dimensional inviscid 
flow and uniform quasi-steady states (1) to (6). In the un- 
choked regime, steady-flow energy, mass, and momentum equa- 
tions were applied across the screen with a screen-drag coefficient 
taken from Adler’s results‘ used in the last. Then application 
of usual shock relations across S$; and Sp and equating of pres- 
sure and velocity in states (2) and (3) permitted iterative solu- 
tion for shock pressure ratios P2,; and Ps. In the choked regime, 
the screens were assumed to be choked isentropic converging- 
diverging nozzles with throat area equal to actual screen open 
area. Then J; = Mscr was readily computed from screen 
blockage, and y = 1.4. Usual shock relations applied across 
S; and Sr then gave P;, from Mscr and Ps. For P22, super- 
sonic state (4) was known from Mscr, P55 and Pa. Shock rela- 
tions applied across Se and Sp with equating of pressure and ve- 
locity in states (2) and (3) then permitted solution for P3, and 

Typical comparisons of theory with experiment are shown for 
the 8-mesh screen in Figs. 2 and 3. The agreement is good con- 
sidering the simplifying assumptions made. 

Attempts to obtain the unchoked screen drag coefficients from 
the measurements met with only limited success. The calcula- 
tion of screen drag from the measured flow quantities was very 
sensitive to experimental inaccuracies and scatter. Further de- 
tails will be found in reference 5. 
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On the Mach Number at the Diffuser Thro 
of an Ejector According to the Hydraulic 
Analogy 


Syogo Matsunaga 


Assistant Professor, University of Osaka Prefecture, Sakai, Osaks 
Japan 
August 8, 1957 


H G. ELROD has studied the Mach Number at the diffyy- 
¢ throat of the ejector in his paper! in which the condition ¢ 
maximizing the induced fluid quantity was used, and when y 
assume the diffuser efficiency equal to 100 per cent, the Mac 
Number at the diffuser throat becomes 7, = 1. 
F Here the author wishes to study the Mach Number at th 
diffuser throat of the ejector by the hydraulic analogy metho 
under the condition of minimizing the absolute pressure of the 
induced fluid 
quantity being suctioned into the mixing chamber. 

The experimental method can be seen in Fig. 1, where 


i.e., in such a condition that there is no fluj 


WS = water supply 

hy = water depth before the nozzle (mm. ) 

hy = water depth at the inlet of the mixing chamber (mm 
hy = water depth at the throat of the diffuser (mm. ) 

hs = water depth at the diffuser exit (mm.) 

R = resistance bar for adjusting the water depth, Hy, 

Q = fluid flow quantity, cu.cm./sec. 

JC = tank for measuring the fluid flow quantity 


In Fig. 2, we can see an example of the experimental results of 
the fluid flow pattern in the ejector type channel under the flow 
condition of minimizing absolute pressure /2 at the inlet of the 
mixing chamber. The disturbances cannot be seen because of 
their small magnitude at the throat of diffuser (a and b in Fig 
2), and the shock waves at the edge of throat are shown in Fig. 2 


Fic. 1. Experimental method. 
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TABLE 1 


hy hy h; hy Q cu.cm./sec. 
A 56 mm. 11 30 25 406 
B 56 8 31 24 400 
B 55 6 30 23 374 

TABLE 2 
pi ata ata 

A 3.9 0.14 0.82 

B 3.3 0.07 0.86 

B 3.4 0.04 0.83 


These flow patterns satisfy the condition of minimizing suction 
absolute pressure—i.e., 2. An example of the experimental 


results is given in Table 1. 

By using the hydraulic analogy relations 

pi/ps = (hi/hs)*?,  po/ps = (he/hs)? 

the absolute pressures p; and p: can be obtained for the delivery 
pressure pj = 1 ata. The Mach Number at the diffuser throat 
(a and b in Fig. 2) can be calculated using the sonic velocity 
mean fluid flow velocity at the throat of diffuser. 

The ratios between the values of hy, he, and h;—i.e., p; and p» 
do not change their values in other experiments with the same 
flow pattern but different magnitudes of 4;, he, and h; in the hy- 


C= gh, (where g is the gravitational acceleration) and the 


draulic analogy. 

The Mach Number ./, at the points a and b in Fig. 2 is shown 
in Table 2. 

The conclusions of the experimental results of hydraulic anal- 
ogy test can be stated as follows: 

To minimize the suction absolute pressure, the Mach Number 
at points a and b must satisfy 17; < 1, and so, there must exit a 
point of 17, = 1 at the inlet into the throat of the diffuser. 
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Drag Due to Lift for Delta Wings at Supersonic 
Speeds? 


Seymour Lampert* 


Research Engineer, Aeronutronic Systems, Inc., 
Glendale, California 
August 8, 1957 


I THIS NOTE the linear theory is compared with some of the 
experimental results obtained in tests! conducted on delta 
wings at the Jet Propulsion Laboratory. 
The drag due to lift according to the linear theory may be 
expressed as a parabolic variation of C, with ACp or 
ACp Cp Comin KC,;? (1) 
For wings with subsonic leading edges (slope of leading edge + 
slope of Mach waves = tan w/tan uw < 1), the parameter K takes 
the form 
Ke = (cot w/27) [E(k) — k/2] (2) 
if the effects of edge suction are included and 


t The results presented here were one phase of the research carried out at 
the Jet Propulsion Laboratory, California Institute of Technology, under 
Joint Sponsorship of the Department of the Army, Ordnance Corps (Under 
Contract No. DA-04-475-Ord 18), and Department of the Air Force 

* Formerly, Research Engineer, Aerodynamic Research Section of the 
Jet Propulsion Laboratory, Pasadena, Calif. 
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Fic. 1. Drag rise coefficient ACp/K vs. Cy, for A,1 wing 
AR = 0.52 and A,2 wing AR = 0.85 with and without edge-suc- 
tion term included at various Mach Numbers. 


> 


K, = (cot w/27) E(k) (3) 
if edge suction effects are not included. E(k) is the complete 
elliptic integral of the second kind whose modulus k = 
[1 — (tan w/tan w)?]'/*. For wings with supersonic leading edges 
tanw/tanyp> 1 

K = 6/4 = V(M? — 1)/4 (4) 
In Figs. 1 through 3 the linear theory is compared with the re- 
sults of tests made for triangular wings of varying aspect ratio 
over the range of Mach Numbers 1.459 to 4.548. The wings 
designed as A, and Ay, are flat sided wings with elliptic leading 
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Fic.2. Drag rise coefficient ACp/K vs. Cr for A,3 wing AR = 
1.43 and A,4 wing AR = 2.40 with and without edge-suction term 
included at various Mach Numbers. 
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Fic. 3. Drag rise coefficient ACp/K vs. Cr for A,.5, A,5, 
and A,o wings AR = 3.62 for various Mach Numbers and Reyn- 
olds Numbers. 


edges (nose radius = 0.022 in.) and sharp leading edges (semi- 
wedge angle in stream direction = 6°), respectively. These 
wings are 2.3 per cent thick based on root chord dimension. 
The wing designated A,» is a diamond profile, max. thickness at 
50 per cent chord and 5 per cent thick. For detailed descrip- 
tions see reference 1. In Figs. 1 and 2 the comparisons are 
made for wings with subsonic edges and the results are presented 
with and without suction effects included. The higher AR wings 
shown in Figs. 2 and 3 had leading edges which theoretically 
could be supersonic over all or part of the Mach Number range 
of the test, hence the proper theoretical value of K is 8/4. 
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On a Reciprocity Condition for Supersonic 
Flutter 


John W. Miles 

Professor of Engineering, University of California, 
Los Angeles, Calif. 

August 8, 1957 


A SSUMING the harmonic time dependence exp (wf), let the 
deflection of a two-dimensional aerodynamic section of 
length / be expanded in the set of functions 

= exp (inrx/l),n = 0, +1, +2,... (1) 
and let p,(x) denote the corresponding set of pressure distribu- 
tions. We seek a reciprocity relation between the generalized 
forces Qm»n and Quin, Where 


1 
Qmn = f, dx 


and 2, is the complex conjugate of zy. 
The pressure distribution on a two-dimensional supersonic 
flow may be expressed in the form! 


(2) 
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where 
g(x) = exp [—1tkM?x/(M? — 1)] Jo[RMx/(M? 1)] 4 


k being the reduced frequency (k = wl/U, U = air speed) and 
the Mach Number. Substituting Eq. (3) in Eq. (2) vields 


= tn(x) dx f. g(x — dé (3 


Now, from Eq. (1), 


zn(x) = (—)"z,(1 — x) 


! 

= (- yr f 
0 

(—)* Sm(&) dé f. g(x — — x) dx (7 
x 

=(-)" f, — x) dx g(x — dE (Te 
! x 


where Eq. (7b) follows from Eq. (7a) by interchanging the order 
of integration, Eq. (7c) from Eq. (7b) by interchanging ¢ and 
1 — x (i.e., introduce the change of variables & = ] — x’ and x = 
1 — £' and then drop the primes), and Eq. (7d) from Eq. (7c) by 
virtue of Eq. (6). Comparing Eq. (7d) to Eq. (5) then yield: 
the reciprocity relation 


whence 


— x) dx g(x — E)Zm(E) dE (7a 


Qmn = (8 


The result (8) also holds if z,(x) = sin (mzx/1), n = 1, 2, 3,.. 
or 3,(x) = cos (mrx/1l),n = 1, 2,3,.... If, on the other hand, 
the sines and cosines are mixed, so that z, = sin (nzx//) but s, = 
cos (mxx/l), the resulting Eq. (8) must be replaced by 


Qmn = (—)™*"*! Qam (9 


The reciprocity relation (8) is of particular interest in panel 
flutter. If the panel is simply supported at x = O and /, the z, = 
sin (wzx/l) form an appropriate set of deflection functions, and 
the coupling coefficients then are nonconservative (Qn, = 
—Qnm) only if m + nis odd. This contradicts Shen’s result*— 
evidently containing an algebraic error in its derivation—that 
Qin = —Qnm for all m and n. The only result actually used in 
Shen’s subsequent analysis was Qi: = —Q, however, so that his 
end results were not affected by the error. 
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An Approximate Solution of the Laminar Heat 
Transfer Along a Heated Flat Plate With an 
Arbitrary Distribution of Surface 
Temperature 


G. Lowe 
School of Engineering, Auckland University College, New Zealand 
August 15, 1957 


V ARIOUS METHODS!~* have been presented for calculating 
laminar heat transfer along a body with nonuniform surface 
temperature. In references 1, 5, and 6, the temperature distribu- 


tion considered is of the type Bx? and the velocity outside the 
boundary layer taken as equal to Cx”. 


In references 3 and 4, 
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ions are applicable to an arbitrary distribution of surface 


the solut 
temperature and stream velocity. Chapman and Rubesin? con- 


sidered uniform flow of a compressible fluid along a flat plate; 
method enabled the heat transfer to be evaluated once the 
nperature distribution along the surface was approximated by 


their 
tel 
a polynominal. 

In this note, a simple approximate method is presented for cal- 


culating the heat transfer along a flat plate directly from the sur- 
face temperature distribution. This method would be useful in 
cases Where the flow is uniform and incompressible. The solu- 
tion is obtained from the well known “heat flow equation” of the 


thermal boundary layer 
C,(d/dx) f, pu(T — Tao) dy = —k(dt/dy)o (1) 


where w is the velocity in the dynamical boundary layer, 7 the 
temperature in the thermal boundary layer. 

In this analysis, the usual assumptions for incompressible flow 
aremade. A solution for Eq. 1 can be obtained by assuming both 
velocity and temperature distributions in the two boundary layers 


to be fourth-degree polynomials 


= 2(¥/5) — 2(y/6)® + (9/6)! 
T/To = 1 + t*(x*){1 — 2(y/51) + 2( (2) 


t*(x*) = t(x*)/T. 
where 6, is the thermal-layer thickness and can differ from the 
dynamical-layer thickness 6, and ¢(x*) the surface temperature 
distribution. The above distributions can be made if the surface 
temperature is not sufficiently large to affect the property values 
and if only positive heat transfer is considered. 
Using Eqs. (2) in Eq. (1), an approximate solution for 


& (= 6,/6) is 


where o is the Prandtl Number. It can be shown that at x* = 0, 
t ~ 0 if the entire plate is heated. Now the local heat-transfer 
coefficient and the local heat transfer are obtained from 


h, = 2k/t, Nu, V Re = ().343/é (4) 


The solution for Eq. 1 was obtained assuming that & did not 
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Fic. 1. Variation of the local heat transfer with ‘‘y’’ for a flat 
plate with surface temperature distribution ¢(x*) = Bx*?, 


READERS’ 


FORUM 


t(x") = 1+ 10x" - 


| 
tix") 


T 
0 40 

| 


| 
—— PRESENTED METHOD _ 
----- CHAP. & RUB.,(REF. 2) 


t 


Fic. 2. Local heat transfer along a heated flat plate with 
approximately constant heat-transfer rate along surface. 


greatly exceed unity. However, for &(x*) = a constant (the 
isothermal plate), & = 0.958 o'/3 = 1,08 for o = 0.70 and gives 
Nu/V Re = 0.318 which is within 8 per cent of Pohlhausen’s 


exact solution. For temperature distribution of the type {(x*) = 


Bx*?, Nu/V Re = 0.3580'/3(2y + 1)'”3 which on comparison 
with the results of reference 5 over the range of 7 = 1 to 4 is within 


8.1 per cent for o = 0.70 and 4.1 per cent for o = 20 (Fig. 1). 

Consider now a heated flat plate with a temperature distribution 
similar to that which would be obtained if there were constant heat 
transfer along the surface (see Fig. 2). (Note that the temperature 
atx* = 1 is six times the leading-edge temperature.) Comparison 
with the present method, Eqs. 3 and 4, with that of Chapman 
and Rubesin’s method shows general agreement for o = 0.72. 
So far, temperature distributions with positive gradients only 
have been considered and, therefore, all the heat transfer is posi- 
tive. For a plate with heat input upstream only, negative heat 
transfer can occur downstream?—a condition in which the present 
method would be unreliable, as temperatures in the boundary 
layer would be greater than at the plate surface. However, for 
temperature distributions with positive gradients only, the pres- 
ent approximate method gives values for the heat transfer com- 
parable with those obtained from exact solutions. While this is 
not a verification of the individual assumption made in the analy- 
sis, it does, however, indicate the usefulness of the final result 
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Electrical Analogy for Transient 
Axisymmetrical Heat Flow 


J.S. Przemieniecki 

Section Leader, Structures Department, Bristol Aircraft Limited, 
Filton, England 

August 15, 1957 


— ELECTRICAL ANALOG for transient heat flow based on the 
fundamental similarity between the flow of heat within a 
homogeneous body and that of charge in a noninductive electric 
circuit has been discussed by Paschkis and Baker.! The exten- 
sion of the method to problems involving axisymmetrical heat 
flow requires only a varying distributed resistance and capaci- 
tance in the one-dimensional analog circuit. Such an analog 
might find useful application in the investigation of temperature 
distributions in circular frames or bulkhead diaphragms in air- 
craft or missile structures subjected to aerodynamic heating.” 
Consider a noninductive hypothetical RC-cable whose dis- 
tributed resistance R’ and distributed capacitance C’ are func- 
tions of the distance x measured from some convenient origin. 
From the fundamental relations for the flow of the electric 
charge, 6Q, in a small element of length, 6x, we have R’éx = 6V/I 
= 6V/(0Q/ot) and C’ix = 6Q/V where V is the voltage and 


J is the current in the cable at point x and time ¢t. Hence, 
= = R’'0Q/dt (1) 
and 
0Q/ox = C’'V (2) 


Differentiating Eq. (1) with respect to x and Eq. (2) with respect 
to ¢, and then eliminating Q we obtain 


(0?V/dx?) — (1/R’) (OR’/dx) (OV/dx) — R'C'(OV/dt) = 0 


(3) 
Now if we make 
—(1/R’) (OR’/dx) = (1/x) (4) 
and 
R’'C’ = (constant) (5) 


then Eq. (3) becomes 
(0?V/ox?) + (1/x) (OV/Ox) — (1/x) (OV/dt) = O (6) 


which is completely analogous to the equation for the conduction 
of heat in problems with axial symmetry.* Namely 


+ (1/r) (OT /dr) — (1/xK) (OT /dt) = (7) 


where 7 is the temperature, x thermal diffusivity and r is meas- 
ured from the center of symmetry. From Eqs. (4) and (5) 


R’ = Ry/x (8) 
and 
C’ = x/KRy (9) 


where Ry is a constant resistance. Thus the necessary require- 
ment for the analog is that the distributed resistance R’ is 
inversely proportional to the distance x, while the distributed 
capacitance C’ is directly proportional to x, where x is measured 
from the center of symmetry. In practice, the hypothetical 
RC-cable is replaced by a “lumped” analog circuit, which 
amounts to replacing distributed resistance and capacitance, as 
given in Eqs. (8) and (9), by a series of varying resistances in 
series and corresponding condensers in parallel. 
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The Properties of Metals Under Rapid Heating 
Conditions 


Elbridge Z. Stowell 
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August 30, 1957 


ia ACCORDANCE with a recently proposed relation for metals, 
an analytical solution has been obtained for the so-called 
“rapid-heating’’ curves. When a metal, kept under a constan 
dead weight stress o, is heated at a constant rate of increase of 
temperature, the solution for the strain « at any absolute tem. 
perature TJ is 


= (0/E) + a-AT + f(AH/RT) >o 


where 
o = stress, ksi 
E = Young's modulus at the temperature 7’, ksi 
a = expansion coefficient at the temperature 7, per °K 
AT = the increase in temperature, °K. 
To = the rate of increase of temperature, °K. per hour 
s = a viscosity constant, per hour per °K. 
AH = a viscosity constant, cal. per mol 
a = aviscosity constant, ksi 
R = the gas constant, 2 cal. per mol per °K. 
f(SH/RT) = {1 — [3/(AH/RT)] + 


[12/(AH/RT)?] — ...} 


The first term gives the contribution due to application of 
stress, the second term gives the contribution due to thermal 
expansion, and the last term gives the contribution due to vis- 
cosity. Since all of the effects of temperature in the final term 
are concentrated in the function f, this function may be tabu- 
lated once and for all in terms of the nondimensional parameter 
(AH/RT). Naturally such a table greatly expedites calculations 
involving rapid heating. 

As an illustration of the application of this solution, Fig. 1 
shows the comparison of the strain-temperature curves for one 
form of Inconel-X, obtained under rapid heating conditions. 
The viscosity constants are s = 10° per hour per °K., AH = 
115,000 cal. per mol, o = 8 ksi. The data were taken 
from reference 2 and apply to Inconel-X sheet which has been 
heat-treated for 1 hour at 1,400°F. and then air cooled. 

The ‘‘yield temperatures’ 7 may be readily obtained from the 
solution given, since for a plastic strain of 0.002 


log f(AH/RTy) = log[7'o/s(AH/R)*] — (0/09) — 6.2 
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© $The right-hand side may be computed from known constants, and 
per hour if log f has been tabulated along with f, the values of AH/RT, 
may be read from the table. The yield temperature 7, is thus 
determined. By way of illustration, Fig. 2 shows the corre- 
sponding comparison of experimental and calculated yield tem- 
peratures for the same material as in Fig. 1, plotted against tem- 
. perature rate. 
| + These solutions are of course only applicable to those metals 
as,” which behave in a sufficiently regular manner so that the vis- 
cosity constants may be determined from measurements of steady 
creep. Once these constants (s, 4H, oo) are determined, the 
other properties of the metal, such as the stress-strain curves, the 
life to rupture, the ultimate strength, the properties when heated 
rapidly, the behavior under conditions of restraint, ete., may be 
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Effect of Externally Generated Vorticity on 
2 Laminar Heat Transfer 


Richard M. Mark 
—— Senior Aerodynamics Engineer, Aerophysics Group, Convair 
Astronautics, San Diego, Calif 


August 30, 1957 


F BLUNT-NOSED BODIES flying at speeds greater than the 
ambient sound speed, vorticity is generated immediately 
behind the curved detached shock wave encompassing the body 
rhe magnitude of this vorticity can become large in the forward 
stagnation region, especially under hypersonic conditions, as 
demonstrated by Lighthill' and can become of the same order 
of magnitude as that generated at the surface of the body. This 
Situation poses a difficult interaction problem between the re- 
spective vorticities generated by the shock wave and the surface. 
However, the usual boundary-layer theory may be applied in a 
modified form to obtain at least a first approximation of the 
effect of externally-generated vorticity if the flow field can be 
divided into two distinct layers (an inviscid rotational layer 
superimposed upon a thin viscous layer) and satisfies certain 
conditions at the dividing surface as pointed out by Ferri and 
‘t for Libby.2 This physical model will be postulated in the present 
irious analysis. 


FORUM 923 


It is the purpose of this note to discuss solely the effect of ex- 
ternal vorticity, such as that generated by a curved shock wave, 
on the laminar heat transfer near the stagnation point of blunt 
bodies of revolution. The following analysis will be simplified 
to bring out the external vorticity effect clearly without sacrificing 
the physical features of the problem 

For flows with small enthalpy differences the density p, co- 
efficient of viscosity u, and the heat conductivity k& are taken to 
be constants. Hence, for a Prandtl Number uC,/k of unity 
(C, is the specific heat at constant pressure), the governing flow 
equations in the viscous layer near the forward stagnation point 
of a blunt body of revolution are 


[O(ux)/Ox] + [d(vx)/Oy] = O (la) 
pu(Ou/Ox) + pv(du/Oy) = —(Op/dx) + (1b) 
pu(Oh/Ox) + prv(dh/Oy) = pw(O7h/dy?) (1c) 


where h is the enthalpy, p = p(x) is the pressure, r = pu(Ou/dy) 
is the shear stress, and w and v are the velocity components in 
the x and y directions (along and normal to the body) respec- 
tively. 

The associated boundary conditions are (a): u = v = O, 
T = Tw, (Op/dx) = (07/dy), and h = h, (constant) at y = 0; 
and (b): u = u,, Ou/Ov = (Ou/Oy), or r = 7, andh = h, (con- 
stant) as y— ©, where the subscripts w and e refer to conditions 
at the wall and the dividing surface respectively. Conditions 
(b) may be obtained by evaluating the results from an inviscid 
flow theory at the wall if the viscous layer is sufficiently thin. 

Since the shear stress varies linearly with x near the stagnation 
point, it can be easily shown that the relation 
h + [(he — hw)/(tw — = Iw + [Che — Ihw)/( te — 

(2) 
satisfies Eqs. (1) and the boundary conditions directly. Hence 
the heat transfer at the wall is proportional to 

[(he — hw)/(tw — [—(Op/Ox)] (3) 

In order to solve for 7,—the only unknown in Eq. (3)—it is 
necessary to transform the governing equations. Combining 
Eqs. (la) and (1b) and using the Crocco variables gives 
— pp u(0/Ox)(x/7r) + wx(Op/dx) (0/du) (1/7) = 

x(0?r/Ou?) (4) 
The associated boundary conditions are 7(07/Ou) = pu(Op/dx) 
atu = Oandr = r,at u = 

Very near the stagnation point, the conditions at the dividing 
surface are written as 
Ue = Vx/a), Te = (Op/Ox) = 
Then a similar solution to Eq (4) exists of the form 

= (pu)? (V/a)*!? xf(n) 
where » = au/Vx, and f is to satisfy the following nondimen- 
sional system 
= (ay? — *) (df/dn) (5a) 
f(df/dn) = —a at n = Oand f = at = a (5b) 
where V is the free-stream velocity, a is the nose radius of the 
body, R is the Reynolds Number (= pVa/u) and the a;,’s are 
dimensionless constants. 

A numerical example will now be presented to illustrate the 
effect of shock-generated vorticity on the heat transfer using the 
inviscid flow results of Lighthill.!. According to Lighthill's re- 
sults, the a;’s are functions of K—the ratio of the density behind 
the shock to that ahead of it. For A = 13.670, a; = 0.418, and 
a, = 10.075. If R = 8.25 X 104, then a,R~'? = 0.035. 

In general, for a; less than unity, system (5) may be solved for 


f—and hence r,—by assuming a power series solution of the form 


f= ; (a,/n!)n", where the coefficients a, may be obtained in 
n=0 

terms of fo = f(m = 0) by substituting the series formally into 


Eq. (5a) and equating coefficients of like powers of » and also 
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using the boundary condition at the wall. The first six coeffi- 


cients are thus: 


a = fo 

a = 

a= —ay*fy 3 

a; = 

ay = —15 7 

a; = 20 5 —105 af, 


A good approximation for f, may be obtained by taking the first 
six terms of the power series as the solution and finding fy that 
will make the resulting relation satisfy the boundary condition 
at 7 = a. In the 
case where the external vorticity is zero or (Ou/Oy), = O, the 
value of fy is 0.855 as computed from the work of Homann.’ 
Thus, with all other factors held constant in Eq. (3), a computa- 
tion of the ratio of the heat transfer in the two cases shows that 
the effect of shock-generated vorticity is to increase the heat 
transfer near the stagnation point by 13 per cent. 

The physical model which was postulated earlier will now be 
If the thickness of the viscous 


In this manner, f) was found to be 0.350. 


justified for the above example. 


Ue 
6 = f (u/7r)du 
0 


then 6 is 0.160 times the shock-detachment distance as computed 


layer is defined by 


One-Dimensional Transient Heat Conduction 
Into a Double-Layer Slab Subjected to a 
Linear Heat Input for a Small Time Interval 


B. Wassermann 

Mathematician, Missiles and Ordnance Systems Department, 
General Electric Co. 

August 12, 1957 


SUMMARY 


a solutions are presented for the transient heat con- 
duction in composite slabs exposed at one surface to a 
triangular heat rate. This type of heating rate may occur, for 
example, during aerodynamic heating. 


SYMBOLS 
b = thickness of first layer, ft. 
| = total plate thickness, ft. 


T = temperature, °F. 


Va 
V/s ky 


ti(x, s) 
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by Lighthill’s theory. Therefore, a thin viscous layer adjacey 
to the surface is distinguished from the rest of the flow field. 7; 
result for fy shows that the shear force at the surface is ten tims 
larger than that at the dividing surface, » = a, and hence ¢} 
dominant shear force is still confined to the layer adjacent to 4 
surface. In the region adjacent to the shock wave the ratio 

the viscous force to the inertia force is of order a2*R™~! since + 

shock detachment distance, x, “4, v, and a» is of order a/K, 

V, V/K, and K respectively. It may thus be concluded that th 
layer is inviscid since a.?R~! is exceedingly small in the aboy 
example. 

The above analysis shows that the effect of shock-generate 
vorticity can become significant in blunt-body heat  transf- 
calculations in the stagnation point region when the paramete 
a2R~'? becomes large, of the order of 0.1 say, and where tj 
Reynolds Number FR in this case is based on the gas properti« 
immediately behind the shock wave, the velocity V ahead 


the shock wave, and the nose radius a. 
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a = fraction of total time, dimensionless 
6 = time, sec. 
Omax = maximum rate of heat input, B.t.u. ft.2/sec 
= heat flux, B.t.u. /ft.?/sec. 
x = distance along x axis, ft. 

= (k2/ki) V ai/a2, dimensionless 
A = (y + 1)/(y — 1), dimensionless 

Subscripts 
1 = initial layer 
2 = secondary layer 
DISCUSSION 


The partial differential equation for transient one-dimensional 


heat flow is 
OT /06 = a(0?7T/dx?) (1 


We assume thermal properties independent of temperature 
and an unheated surface insulated; at the interface, there is m 
resistance to heat flow. The initial temperature is, of course 
0 at 0 = O (see Fig. 1). 

Solving for the temperature, one obtains, 


J (y+1) cosh V s/ara: [Va = b) + Va: x) — (7 — 1) cosh WV a, (L — b) - 


(y + 1) sinh V (L—b) + Vaz | + (y — 1) sinh (L — b) - V as 


tidt ale 
s) = a 
sky 


where / is the transform of J and s the transform of 6. 

In order to invert the Laplace transform, the relative sensi- 
tivities of the two layers are restricted to three cases, respec- 
tively,! 

= = = 
ki/Vas = ko/V ar, a2 > ki/V a, and 
ki/V a, > Va 


The temperatures in the first and second layers for Case I are 


_2 cosh V's/az (L x) 


(y + 1) sinh V s/aya: (Va, (L — b) + V a2 | + —1)sinh Ys (L — - Va 


Ti(x, 6) = 473 erfe 
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2(n + l)s V a2 xl (4) 
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where = Va 6) + V a2 b 


/ 
V a 


2nz + a2Xx 
n(a, 0) = 52 erfe 


where 2 = Va (L — b) + V/ a: b. For Case II, 
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aR n=0 
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where 


The heat flux at any time and position may be obtained by 
first differentiating Eqs. (2) and (3), and then inverting the 


(6) 


(2n + 1)b+(b —x) 


(2n + 1)b — (h — x) 


(7) 


The solution for the composite semi-infinite slab, becomes, for 


(9) 


be 

Fic. 1. Right: One-dimensional heat flow in a double-layer 
slab. Left: Heat input 
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ma 1 2/9 
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7 (10) 
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Eq. (10) was derived? by considering a semi-infinite composite 


slab. For Case III, 


ah; n=0 
(2n +1)b —(L (2n+1)b+(L —x) 
+ Perfc . 
2V a0 2V a6 { 
(9 ad) erfe : 
(1 — a) n=0 2V a(@ — af) 
2n + 1)b +(L — x)) 
erfe €11) 


/ 
2V «(6 — 
Since the secondary layer acts as an insulator the heat transfer 
across the interface can be considered almost nil 
It must be remembered that if 6/06 < a, then (0/6) — a = 0. 
The Laplace transform of Fig. 1. (Jeft) is 


= Qmax(1 — a — e~@5)/(1 — a)ads? (12) 


Transforms of other types® of linear heat rates can also be sub- 
stituted for L}q{ before applying the inverse transformation, 


APPENDIX 


In these formulas the repeated integrals of the complementary 
error function is* 


i" erfex = f erfe dyn 
x 


Then integrating by parts we have 


ierfex = (lV — xerfex 
The general! recurrence tormula is 
2ni" erfe x = erfe x — 2xi"~' erfe x 
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Manuscripts: The original typewritten copy of the paper is 
desired, double or triple spaced on one side of white paper sheets, 
consecutively numbered. There should be wide margins to allow 
for the marking of directions to the printer. Correcting, chang- 
ing, or adding to papers after they are in type is costly. It is, 
therefore, imperative that papers submitted be in final form. 
Typographical errors may be corrected on proofs, but if authors 
wish to add material, they may do so at their own expense. In 
mailing, drawings may be rolled, but manuscripts should be sent 
flat. Send by first-class mail (register if you wish for your own 
protection) to the Editorial Office, Institute of the Aeronautical 
Sciences, 2 East 64th Street, New York 21, N.Y. All manuscripts 
will be examined by the Editorial Committee and by the Editor. 
Authors will be advised as promptly as possible whether the paper 
is acceptable for publication. 


TitLEs: The title of the paper should be brief. The name 
and initials of the author should be written as he prefers. The 
use of the full name of an author is advocated because of the fre- 
quent duplication of initials and surnames which sometimes makes 
it difficult to establish the identity of the author. This is par- 
ticularly important for large annual indexing and abstracting 
services. All titles and degrees or honors are omitted. The name 
of the organization with which the author is associated should 
be placed after his name on a separate line. The date on which 
the paper is received will be inserted by the Editor. The author’s 
title or position should be indicated in a footnote. 


SumMMARIES OR ABsTRACTS: An abstract to be printed at the 
beginning should accompany each article. It should be adequate 
as an index and as a summary. It should contain a statement of 
major conclusions reached, since summaries in many cases con- 
stitute the only source of information used in compiling scien- 
tific reference indexes. Abstracts printed in other journals, espe- 
cially foreign, in most cases, consist of summaries from printed 
papers. The summary should explain as adequately as possible 
the major conclusions to a nonspecialist in the subject and should 
contain from 100 to 300 words, depending on the length of the 


paper. 


SuBHEADINGS: Subheadings should be inserted by the author 
at frequent intervals. The work of editorial preparation will be 
simplified by the author’s provision of many subheadings. 


MatTTER UsuaL_ty DELETED: Photographs or illustrations of 
little technical interest and not showing advances in general 
practice. Too detailed tabular matter (general results of such 
tables may be included in the text). Lengthy descriptions of 
materials or processes or of preliminary experiments or theories 
that preceded final results; salient features only are of interest. 


REFERENCES: References should be numbered consecutively 
and grouped together at the end of the manuscript, with only 
the corresponding number being mentioned im- the text. The 
arrangement should be as follows. For books: ! Durand, W. F., 
Aerodynamic Theory, 1st Ed., Vol. 1, p. 23; Julius Springer, 
Berlin, 1934. For magazines: ' England, C. R., Crawford, 
A. B., and Mumford, W. W., Some Results of a Study of Ultra- 
Short-Wave Transmission Phenomenon, Proc. I.R.E., Vol. 20, 
No. 12, pp. 481-482, March, 1933. Please give author, title, 
edition, volume, number, page, publisher, and date of publication 
as indicated. Omission of one required fact causes much extra 
editorial work and possible inaccuracies. 


ILLUSTRATIONS: Illustrations should accompany manuscripts, 
and each should always be referred to in the text by number. 
Drawings or graphs should not be larger than 12 X 16 inches, 
and must be made with jet black India ink on white paper or 
tracing cloth, the latter being preferred. Do not use typewriter 
for lettering. The smallest lettering on 8 X 10 inch figures should 
be no less than '/, inch high. Cross-section paper (white with 
black lines) may be used, but it should not have more than 4 lines 
per inch. If finer ruled paper is used, the major division lines 
should be drawn in with black ink, omitting the finer divisions. 
In the case of finely ruled paper, only blue-lined paper can be 
accepted. Tracing paper and blueprints are not acceptable. 
Lettering and all markings must be large enough to be readable 
after reduction to single-column width (3°/,.in.). Mail rolled or 
flat; never fold. Drawings that cannot be reproduced (including 
pencil drawings) will be returned to the author for redrawing, 
thus delaying publication of the paper. Photographs should be 
distinct and show clear black and white contrasts. They must 
be on glossy white paper. Avoid round and oval photographs. 


CAPTIONS AND LEGENDS: Legends or captions must accompany 
each drawing or photograph submitted. If written on the draw- 
ing or photograph, they should be placed below and well outside 
the part to be reproduced. Each table should have a caption 
such as Table 1, Table 2, Table 3, etc. Captions should be com- 
plete in themselves so as to make the data intelligible to the 
reader without reference to the text. A duplicate list of captions 
for figures should be included as the last page of the manuscript. 
Use “Fig. 1” (not Figure 1), ‘‘Figs. 3 and 4,”’ etc., in both the text 
and the numbering of illustrations. In the text, ‘Eq. (1)” or 
“Eqs. (1) and (2)” is used, not “Equation (1). In captions, 
legends, and in table headings, write all words in full; do not 
abbreviate, except for “‘Fig.’’ and “‘Eq.” 


MATHEMATICAL WorRK: Only the simplest formulas should be 
typewritten; all others should be carefully written in pen and 
ink, the writing to be large enough so that it can be marked 
for the printer. Considerable space for marking should be al- 
lowed above and below all equations. All complicated equa- 
tions should be repeated on separate sheets with plenty of 
space left for marking. The solidus should be used for sim- 
ple fractions appearing within the text. Make all expressions 
clear to the typesetter. Greek letters used in formulas should 
be clearly designated by name in the margin of the manu- 
script. The difference between capital and lower-case letters 
should be clearly distinguished and care taken to avoid confu- 
sion between zero (0) and the letter (0), between the numeral 
(one) and the letter (ell) and the prime (’), between alpha and 
a, kappa and k, u and mu, v and nu, n and eta. All sub- 
scripts and exponents should be clearly distinguished. Avoid 
complicated exponents and subscripts. Dots and bars over 
letters or mathematical expressions should also be avoided. 
When it is necessary to repeat a complicated expression, it should 
be represented by some convenient symbol. 


SYMBOLS AND ABBREVIATIONS: The symbols recommended in 
the American Standards Association Standard “Letter Symbols 
for Aeronautical Sciences,’’” ASA Y10-7—1954, should be used 
wherever practicable. All symbols should be clearly written and 
carefully checked. Standard abbreviations should be used, and 
it should be noted that most abbreviations are lower case, such 
as m.p.h., b.m.e.p., ih.p., b-hp., hp., etc. 
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